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Abstract 

We  show  that  Jacobi's  method  (with  a  modified  stopping  criterion)  computes  small 
eigenvalues  of  symmetric  positive  definite  matrices  with  uniformly  higher  relative  accu- 
racy than  QR,  divide  and  conquer,  traditional  bisection,  or  any  algorithm  which  first 
involves  tridiagonalizing  the  matrix.  If  fact,  modulo  an  assumption  based  on  extensive 
numerical  tests,  we  show  that  Jacobi's  method  is  optimally  accurate  in  the  following 
sense:  if  the  matrix  entries  have  small  relative  errors,  so  that  the  eigenvalues  have  cer- 
tain intrinsic  uncertainties,  Jacobi  will  compute  them  with  nearly  this  accuracy.  In 
other  words,  as  long  as  the  initial  matrix  has  small  relative  errors  in  each  component, 
even  using  infinite  precision  will  not  improve  on  Jacobi  (modulo  factors  of  dimensional- 
ity). We  do  prove  bisection  is  optimal  in  this  sense.  We  also  show  the  eigenvectors  are 
computed  more  accurately  by  Jacobi  (and  inverse  iteration)  than  previously  thought 
possible.  We  prove  similar  results  for  using  one-sided  Jacobi  for  the  singular  value  de- 
composition of  a  general  matrix.  We  present  a  version  of  Jacobi  with  the  property  that 
the  more  its  accuracy  exceeds  that  of  QR,  the  faster  it  converges. 
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1      Introduction 

Jacobi"s  method  and  QR  iteration  are  two  of  the  most  common  algorithms  for  solving 
eigenvalue  and  singular  value  problems.  Both  are  backward  stable,  and  so  compute  all 
eigenvalues  and  singular  values  with  an  absolute  error  bound  equal  to  p(n)r  H/Z^Ht.  where 
p{n)  is  a  slowly  growing  function  of  the  dimension  n  of  the  matrix  H ,  s  is  the  machine 
precision,  and  \\H\\2  is  the  two  norm  of  the  matrix.  Thus,  large  eigenvalues  and  singular 
values  (those  near  ||i^||2)  are  computed  with  high  relative  accuracy,  but  tiny  ones  may  not 
have  any  relative  accuracy  at  all.  Indeed,  it  is  easy  to  find  symmetric  positive  definite 
tridiagonal  matrices  where  QR  returns  negative  eigenvalues.  This  error  analysis  does  not 
distinguish  Jacobi  and  QR,  and  so  one  might  expect  Jacobi  to  compute  tiny  values  with  as 
little  relative  accuracy  as  QR. 

In  this  paper  we  will  show  that  Jacobi  (with  a  modified  stopping  criterion)  computes 
eigenvalues  of  positive  definite  symmetric  matrices  and  singular  values  of  general  matrices 
with  a  uniformly  better  relative  error  bound  than  QR,  or  any  other  method  which  initially 
tridiagonalizes  (or  bidiagonalizes)  the  matrix.  This  includes  divide  and  conquer  algorithms, 
traditional  bisection  and  Rayleigh  quotient  iteration,  etc.  We  will  also  show  Jacobi  computes 
eigenvectors  and  singular  vectors  with  better  error  bounds  than  QR. 

In  fact,  for  the  symmetric  positive  definite  eigenproblem.  we  will  show  that  Jacobi  is 
optimally  accurate  in  the  following  sense.  Suppose  the  initial  matrix  entries  have  small  rel- 
ative uncertainties,  perhaps  from  prior  computations.  The  eigenvalues  will  then  themselves 
have  inherent  uncertainties,  independent  of  which  algorithm  is  used  to  compute  them.  We 
will  show  that  the  eigenvalues  computed  by  Jacobi  have  error  bounds  which  are  nearly  as 
small  as  these  inherent  uncertainties.  In  other  words,  as  long  as  the  initial  data  is  slightly 
uncertain,  even  using  infinite  precision  cannot  improve  on  Jacobi  (modulo  factors  of  n).  For 
the  singular  value  decomposition,  we  can  prove  a  similar  but  necessarily  somewhat  weaker 
result. 

These  results  depend  on  new  perturbation  theorems  for  eigenvalues  and  eigenvectors  (or 
singular  values  and  singular  vectors)  as  well  as  a  new  error  analysis  of  Jacobi,  all  of  which 
are  stronger  than  their  classical  counterparts.  They  also  depend  on  an  empirical  observa- 
tion for  which  we  have  overwhelming  numerical  evidence  but  somewhat  weaker  theoretical 
understanding. 

Now  we  discuss  our  error  bounds  in  more  detail.  First  we  consider  the  standard  error 
bounds.  Let  H  he  a,  positive  definite  symmetric  matrix,  and  SH  a  small  perturbation  of  H 
in  the  sense  that  [SHtj/Hi^]  <  rj/n  for  all  i  and  j.  Then  ||<5/r||2  <  T/IJ/fUj.  Let  A,  and  X[ 
be  the  i-th  eigenvalues  of  H  and  H  +  6H ,  respectively  (numbered  so  that  Ai  <  •  •  •  <  A„). 
Then  the  standard  perturbation  theory  [12]  says 


A.        -       A. 


where  k(H)  =  \\H\\2  •  \\H    '||2  is  the  condition  number  of  H.  We  wiU  prove  the  following 

1  /2 

stronger  result:  Write  H  =  DAD  where  D  =  diag  (.ff,^    )  and  Au  =  1.   By  a  theorem  of 
Van  der  Sluis  [13,6]  k{A)  is  less  than  n  times  nun^K{Dn D),  i.e.  it  nearly  minimizes  the 


condition  number  of  H  over  all  possible  diagonal  scalings.  Then  the  following  is  true: 

\^l^<r,K(A)  (1.2) 

i.e.  the  error  bound  rjK{H)  is  replaced  by  t?k(/1).  Clearly,  it  is  possible  that  k(.-1)  <  k.{H) 
(and  it  is  always  true  that  k{A)  <  nK(H)).  so  the  new  bound  is  always  at  least  about  as 
good  and  can  be  much  better  than  the  old  bound. 

In  the  case  of  the  singular  values  of  a  general  matrix  G,  we  similarly  replace  the  conven- 
tional relative  error  bound  riK{G)  with  t]k.(B),  where  G  =  BD.  D  chosen  diagonal  so  the 
columns  of  B  have  unit  two-norm.  This  implies  k(B)  <  n^^^  min^K(GD).  and  as  before  it 
is  possible  that  k{B)  <^  k{G). 

We  will  weaken  the  assumption  of  small  componentwise  relative  error  | <5 //■,_, //r,j|  <  rjjn 
to  include  the  more  general  perturbation  \6Htj\/iH,iHjj)^^'^  <  T]/n.  These  more  general 
perturbations  include  the  rounding  errors  introduced  by  Jacobi,  which  will  imply  that  Jacobi 
computes  eigenvalues  with  the  error  bound  0{sk(A)).  In  contrast.  QR  (or  any  algorithm 
which  first  tridiagonalizes  the  matrix)  only  computes  them  with  error  bound  0(£k{H)). 

Now  we  consider  the  eigenvectors  and  singular  vectors.  Again  let  77  be  a  positive 
definite  symmetric  matrix  with  eigenvalues  A,  and  unit  eigenvectors  v,.  Let  6H  he  k  small 
componentwise  relative  perturbation  as  before,  and  let  \[  and  v[  be  the  eigenvalues  and 
eigenvectors  oi  H  +  6H .  Then  the  standard  perturbation  theory  [12]  says 

\\v,-v[\\<--^ +  0{v^)  (1.3) 

absgapx, 

where  the  absolute  gap  for  eigenvalues  is  defined  as 

.     |A, -AJ 
absgapx,  =  nun  —  (1-4) 

j#>      ||77||2 

We  will  prove  a  generally  stronger  result  which  replaces  this  bound  with 

where  the  relative-  gap  for  eigenvalues  is  defined  as 

The  point  is  that  if  H  has  two  or  more  tiny  eigenvalues,  their  absolute  gaps  are  necessarily 
small,  but  their  relative  gaps  may  be  large,  so  that  the  corresponding  eigenvectors  are  really 
well-conditioned.  We  will  prove  an  analogous  perturbation  theorem  for  singular  vectors 
of  general  matrices.  We  will  also  prove  a  perturbation  theorem  which  shows  that  even 
tiny  components  of  eigenvectors  and  singular  vectors  may  be  well-conditioned.  Again,  we 
show  Jacobi  is  capable  of  computing  the  eigenvectors  and  singular  vectors  to  their  inherent 
accuracies,  but  QR  is  not. 


To  illustrate,  consider  the  symmetric  positive  definite  matrix  H  =  DAD  where 


H  = 


10^° 

10-9 

1019  ■ 

1 

.1 

.1 

10'^ 

lO^o 

109 

.    A  = 

.1 

1 

.1 

10'9 

10^ 

1 

.1 

.1 

1 

and       Z)  =  diag(10'",  10'"a) 


Here  k{H)  x  10''°  and  k{A)  s;  1.33.  Thus,  rj  relative  perturbations  in  the  matrix  entries 
only  cause  -iq  relative  pert urbaf-ions  in  the  eigenvalsies  accordiug  to  the  nev;  theorem,  and 
3  ■  10'*°  •  rj  relative  \>erturbations  according  to  the  conventional  theorem.  Also,  the  absolute 
gaps  for  the  eigenvalues  of  H  are  absgapx^  2.3  ~  10~^°,  10"'^°,  1,  whereas  the  relative  gaps 
relgapx-^  2,3  are  all  approximately  10^°.  Thus  the  new  theory  predicts  errors  in  v-i  and  i'2  of 
norm  2-10~i°r7,  whereas  the  old  theory  predicts  errors  of  lO^^rj.  Jacobi  will  attain  these  new- 
error  bounds,  but  in  general  QR  will  not.  For  this  example,  QR  computes  two  out  of  the 
three  eigenvalues  as  negative,  whereas  H  is  positive  definite.  In  contrast,  Jacobi  computes 
tdl  the  eigenvalues  to  nearly  full  machine  precision.  In  fact  for  this  example  we  will  show 
Jacobi  computes  all  components  of  aJl  eigenvectors  to  nearly  full  relative  accuracy,  even 
though  they  vary  by  16  orders  of  magnitude;  again  QR  will  not  even  get  the  signs  of  many 
small  components  correct. 

The  empiriccil  observation  on  which  the  ultimate  accuracy  of  Jacobi  depends  is  as  follows: 
Let  Ho  =  DqAqDo  be  the  original  matrix,  and  Hm  =  DmAmDm  where  Hm  is  obtained  from 
Hm-i  by  applying  a  single  Jacobi  rotation,  Dm  is  diagonal,  and  Am  has  unit  diagonal.  The 
desired  error  bound  is  proportioned  to  k(/1o),  but  Jacobi's  error  bound  is  proportional  to 
max„,K(Am)-  So  for  Jacobi  to  attciin  the  claimed  accuracy,  majCm  «(>!„)/«( /lo)  must  be 
modest  in  size.  In  over  900000  numerical  experiments,  its  maximum  value  was  less  than 
1.82.  Our  theoretical  understanding  of  this  excellent  behavior  is  incomplete  and  providing 
it  remains  an  open  problem. 

We  will  also  show  bisection  and  inverse  iteration  (with  appropriate  pivoting,  and  applied 
to  the  original  positive  definite  symmetric  matrix)  are  capable  of  attaining  the  same  error 
bounds  as  Jacobi.  Of  course  bisection  and  inverse  iteration  on  a  dense  matrix  are  not 
competitive  in  speed  with  Jacobi,  unless  only  one  or  a  few  eigenvalues  are  desired  and  good 
starting  guesses  are  available. 

This  work  is  an  extension  of  work  in  [2],  where  analogous  results  were  proven  for  matrices 
which  are  called  scaled  diagonally  dominant  (s.d.d.).  The  positive  definite  matrix  H  = 
DAD  is  s.d.d.  if  rA  is  diagonally  dominant  in  the  conventional  sense.  This  work  replaces 
the  assumption  that  A  is  diagonally  dominant  with  mere  positive  definiteness,  extending 
the  results  of  [2]  to  all  positive  definite  symmetric  matrices,  as  well  as  to  the  singular  value 
decomposition  of  general  matrices. 

This  work  does  not  contradict  the  results  of  [7,2]  where  it  was  shown  how  a  variation 
of  QR  could  compute  the  singiilar  values  of  a  bidiagonal  matrix  or  the  eigenvalues  of  a 
symmetric  positive  definite  tridiagonal  matrix  with  high  relative  accuracy.  This  is  because 
reducing  a  dense  matrix  to  bidiagoncd  or  tridiagonal  form  can  cause  large  relative  errors  in 
its  singular  values  or  eigenvalues  independent  of  the  accuracy  of  the  subsequent  processing. 
In  contrast,  the  results  in  this  paper  are  for  dense  matrices. 

We  will  also  present  an  accelerated  version  of  Jacobi  for  the  symmetric  positive  definite 
eigenproblem  with  an  attractive  speedup  property:  The  more  its  accuracy  exceeds  that 
attainable  by  QR  or  other  traditional  methods,  the  faster  it  converges.  See  also  [16]  where 


earlier  references  for  Jacobi  methods  on  positive  definite  matrices  as  well  as  for  one  sided 
methods  can  be  found. 

The  rest  of  this  paper  is  organized  as  follows.  Section  2  presents  the  new  perturbations 
theorems.  Section  3  discusses  two-sided  Jacobi  for  the  symmetric  positive  definite  eigen- 
problem..  Section  4  discusses  one-sided  Jacobi  for  the  singular  value  decomposition.  It  also 
presents  the  accelerated  version  of  Jacobi  just  mentioned.  Section  5  discusses  bisection  and 
inverse  iteration.  Section  6  discusses  bounds  on  maxm 'c(.4m  )/«(-4o)-  Section  7  contains 
numerical  experiments.  Section  8  presents  our  conclusions  and  discusses  open  problems. 


2      Perturbation  Theory 

In  this  section  we  prove  new  perturbation  theorems  for  eigenvalues  and  eigenvectors  of 
symmetric  positive  definite  matrices,  and  for  singular  values  and  singular  vectors  of  general 
matrices.  In  the  first  subsection  we  consider  eigendecompositions  of  symmetric  positive 
definite  matrices.  In  the  second  subsection  we  discuss  the  optimaJity  of  these  bounds.  In 
the  third  subsection  we  consider  the  singular  value  decomposition  of  general  matrices.  In 
the  fourth  subsection,  we  discuss  the  optimality  of  this  second  set  of  bounds. 

2.1      Symmetric  Positive  Definite  Matrices 

The  next  two  lemmas  also  appeared  in  [2]: 

Lemma  2.1  Let  H  and  K  he  symmetric  matrices  with  K  positive  definite.  Let  the  pencil 
H  -  XK  have  eigenvalues  A,.  Let  SH  and  6Ii  be  symmetric  perturbations  and  let  \[  be  the 
eigenvalues  of  (H  +  SH)  -  A(  A'  +  SK).  Suppose  that 

\x^6Hx\  <  riH  ■  \x'^nx\      and        \x'^6Kx\  <  t^k  ■  l-t^A'j| 

for  all  vectors  x  and  some  r/^j  <  1  and  rjj^  <  1.    Then  either  A,  =  A^  =  0  or 

1  -  riH  ^  K  ^  1  +  VH 

for  all  i. 

Proof.     The  condition  on  6K  implies  that  K  +  SK  is  positive  definite  too,  so  the 

perturbed  pencil  is  definite.  The  condition  on  SH  implies  that  it  and  H  have  the  same  null 

space,  implying  that  A,  =  A(  =  0  if  one  of  them  equals  zero.    Now  we  consider  the  case 

A,  >  0;  for  the  other  eigenvalues  consider  -H  and  -SH .    The  Courant-Fischer  minimax 

theorem  [12]  expresses  A,  as 

x^Hx 
A,  =  rmn  max    - 

S"  xeS'  "^^^^ 
where  the  minimum  is  over  all  i-dimensional  subspaces  S'.  Let  the  spaces  Sq  and  Sj  satisfy 

•    ,                 x^Hx          ^       ^,               x^iH  +  SH)x 
A,-  =  max    -  ,.        and       A,  =  max  -^=- 

x€s;^^Ax  •    ^^s\^^U<  +  6K)x 

Then 

x^(H  +  6H)x^  x^iH  +  SH)x  x^Kx  x^Hx       l  +  m. 

A   =  min  max  -r-— tjt—  <  max =r— •  -^jr—- tttt—  •     x  r.     <  ^'^t 

S'  xeS"  ^  (^  +  <5^)^      leS;        ^^^  x^(K  +  6K)x    x^Kx  -  I-tj^ 

and  similarly 

x^Hx  ^  x'^Hx  x^(K  +  SK)x    x^jH  +  SH)x        I  +  tjk   „ 

A.  -  ininm|c  ^^^.^  <  m^K^  -^-^——^  .         ^^^^         .  _______  <  ^y— ;^A. 

completing  the  proof.       I 


Lemma  2.2  Let  H  =  A]^.4//A//  and  .A//  be  symmetric  matrices.  H  and  A//  need  not  have 
the  same  dimensions,  and  Ah  may  be  an  arbitrary  full  rank  conforming  matrix.  Similarly. 
let  K  =  A^-rlA'AA'  and  A^  be  symmetric  positive  definite  matrices,  where  K  and  A^- 
need  not  have  the  same  dimensions  and  Ar-  may  be  an  arbitrary  full  rank  conforming 
matrix.  Let  SH  -  AJ^SAhAh  be  a  perturbation  of  H  such  that  \x'^SAhx\  <  r]H\x^--iH^\ 
for  all  X  where  rjfj  <  l-  Similarly,  let  (SA'  =  A]^-6AkAk  be  a  perturbation  of  K  such  that 
\x^6Akx\  <  r)i^-\x^AKx\  for  all  x  where  tjk  <  1.  Let  A,  be  the  i-th  eigenvalue  of  H  -  XK 
and  \[  the  i-th  eigenvalue  of  {H  +  ill)-  \[K  +  <5A').    Then  either  A,  =  A',  =  0  or 


1  +  '7ff  ~  A,  -  1 


riK 


Proof.    Note  that  for  all  vectors  x 


\x'^6Hx\  =  |x^aJ(5.-1hAhx|  <  7?„|x^A^AhA//i|  =  r)H\x'^Hx\ 


and 


\x'^U<x\  =  \x^a16AkAkx\  <  tjkIx^aI^AkAkxI  =  nK\x'^Kx\ 
Now  apply  Lemma  2.L       I 


Theorem  2.3  Let  H  =  DAD  be  a  symmetric  positive  definite  matrix,  and  D  =  diag  (/T,, 


1/2n 


SO  .4,1  =  L    Let  SH  =  D&AD  be  a  perturbation  such  that  \\SA\\2  =  ^  <  -^inm(^) 
the  i-th  eigenvalue  of  H  and  X[  be  the  i-th  eigenvalue  of  H  -\-  SH.    Then 


A,  -  a; 


A, 


< 


<  n(^)-  V 


Aminl  ^  ) 

In  particular,  if\SH,j/nij\  <  rj/n,  then  \\SA\\2  <  rj  and  the  bound  (2.4)  applies. 
Proof.    Note  that  for  all  nonzero  vectors  x 


Let  A,  be 


(2.4) 


:T6H: 


t^Hi 


x'^A'^SAAx 


xTA^AAx 


y'^SAy 


y'^Ay 


< 


.(A) 


Lemma  2.2  yields  the  desired  bound,  using  K  =  I  and  ^A'  =  0.  It  remains  to  prove  that 
\6n,j/ntj\  <  ri/n  implies  ||^A||2  <  V-  But  A„  =  1  and  A  positive  definite  imply  that  no 
entry  of  A  is  larger  than  1  in  absolute  value.  (Note  that  this  means  k(A)  is  at  most  n  times 
larger  than  1/Ainin(^).)  Therefore  \6Aij\  =  \6n,j/H,j  ■  A,j\  <  ri/n  and  so  ||^>1||2  <  ^  as 
desired.       I 

Proposition  2.13  in  the  next  subsection  wiU  show  that  the  bound  of  Theorem  2.3  is 
nearly  attained  for  at  least  one  eigenvalue.  However,  other  eigenvcdues  may  be  much  less 
sensitive  than  this  most  sensitive  one.  The  next  proposition  provides  individual  eigenvalue 
bounds  which  may  be  much  tighter: 


Proposition  2.5   Let  H  =  DAD  be  as  in  Theorem  2.3,  with  eigenvalues  A,  and  unit  eigen- 
vectors I',.    Let  H  +  6H  =  D(A  +  bA)D  have  eigenvalues  AJ.    Let  \\SA\\2  =  rj  <^  Anun(-A)- 

Then  the  bound 

I  \         \'i        ^11  r),..i 

^^+0(v^)  (2.6) 


|A.  -  A;|  ^  v\\Dv,\^^ 
A.        -        A. 


is  attainable  by  the  diagonal  perturbation  6Ajj  =  r/sign(t;,(j)). 


Proof.  Bound  (2.6)  is  derived  from  the  standard  first  order  perturbation  theory  which 
says  XAH  +  ^H)  =  X,{H)  +  vJSH v,  +  0(\\6H\\l).  and  snhst\tuUng\vf SH v,\  =  [vJObADv,]  < 
\\Di\\\l\\SA\\2-  The  inequality  \vj D6ADi\\  <  \\Dvi\\l\\6A\\2  is  clearly  attained  for  the  diag- 
onal choice  of  dA  in  the  statement  of  the  proposition.      I 

We  stated  Lemmas  2.1  and  2.2  separately  in  order  to  emphasize  their  generality.  For 
example,  suppose  H  and  A'  are  a  finite  element  matrices.  A;/  and  A^'  the  fixed  assembly 
matrices  of  O's  and  I's,  and  Ah  and  .4a-  the  block  diagonal  matrices  of  individual  elements. 
If  6Afi  and  6Af^  are  also  block  diagonal,  perturbmg  each  separate  element  in  the  sense 
of  x'^SAh^/^^ Ah^  and  x^6Akx/x^ A^-x  being  small  for  all  nonzero  x.  t'len  the  relative 
perturbations  of  the  eigenvalues  of  the  assembled  matrix  H  will  also  be  small.  In  particular, 
consider  the  matrix  arising  from  modeling  a  series  of  masses  mi, . . . ,  m„  on  a  line  connected 
by  simple  linear  springs  with  spring  constants  kQ,...,kn  (the  ends  of  the  extreme  springs 
are  fixed).  The  natural  frequencies  of  vibration  of  this  system  are  the  square  roots  of  the 
eigenvalues  of  the  pencil  M  —  XK  where  M  =  diag  (mi, ... ,  m„)  and  A'  is  tridiagonal  with 
diagonal  /cq  +  ki,  ki  +  k2,  ■  ■  ■ .  kn-i  +  kn  and  offdiagonal  -fci, . . . ,  -/:„-!•  Lemma  2.2  implies 
that  an  rj  relative  perturbations  in  any  spring  constant  or  mass  makes  at  most  rj  relative 
changes  in  the  eigenvalues  of  M  —  XK .  Unfortunately,  this  property  does  not  extend  to  the 
matrix  assembled  in  floating  point,  since  if  A:,_i  <C  A:,  >  fc,+i  so  that  A;,  +  fc,±i  rounds  to  A:,, 
the  computed  A'  will  be  singular  instead  of  positive  definite,  meaning  that  some  eigenvalues 
are  completely  changed.  ' 

One  may  also  prove  a  version  of  Lemma  2.1  in  an  infinite  dimensional  setting  [10,  section 
VI.3]. 

Now  we  turn  to  eigenvectors.  A  weaker  version  of  the  following  theorem  also  appeared 
in  [2]: 

Theorem  2.7  Let  H  =  DAD  be  as  in  Theorem  2.3.  Define  n{t)  =  D(A-\-  tE)D,  where 
E  is  any  matrix  with  unit  two-norm.  Let  A,(e)  be  the  i-th  eigenvalue  of  H{e),  and  assume 
A,(0)  is  simple  so  that  the  corresponding  unit  eigenvector  v,(€)  is  well  defined  for  sufficiently 
small  e.   Then 

||,,„  _  „,(0,||,  <         '" -')■;''        +  0,e',  <  '"-''■"-'^"  +  OU^) 
>^min(A)  ■  relgapx,  relgapx, 

Proof.    Let  vjt(O)  be  abbreviated  by  Vk.  From  [9]  we  have 

t;.(e)  =  v,+€Y^  -\ — i  .  Vk  +  0(€2) 


Let  yk  =  Dvk,  so  that 


v,{e)  =  v.  +  e  ^  |i^  •  Vk  +  0{e^)  (2.8) 

i^  A.  -  Afc 


The  pair  (A,,  y,)  is  an  eigenpair  of  the  pencil  A  -  XD   ^ .  Thus 

>^k  =  XkyjD'^yk  =  VkAyk  >  X^niA)\\yk\\l 


and  so  \\yk\\2  <  ( Afc/A^„(  .4))'/2    Letting  z^  =  yk/Wykh  lets  us  write 

(A, -A,)/(AfcAji/ 


^M  =  v,+eZ,,h::lf~:,,^.-^^^oi^) 


where  |^.i|  =  Ht/fclbilyJb/i AfcA,)^/'  <  1/A^n(.4).  Taking  norms  yields  the  result.       I 

Proposition  2.14  in  the  next  subsection  will  show  that  the  bound  in  Theorem  2.7  is 

nearly  attainable  for  all  v,. 

As  in  Corollary  3  in  [2],  it  is  possible  to  derive  a  nonasymptotic  result  from  Theorem 

2.7: 

Corollary  2.9   Let  H  =  DAD  be  as  in  Theorem  2.3.  Suppose  6  =  ||<5.4||2/Aniin(.A)  satisfies 

.       1  ,       3  ■  2-1/2 -^ 

6  <  -     and       ■ —  <  relgap\, 

4  1  —  0 

Let  u,   be  the  i-th  unit  eigenvector  of  H  =  DAD  and  v[  be  the  i-th  unit  eigenvector  of 
H'  =  D{A^  8.A)D.   Then 

[n-lY'^b 


:i  -  Ab){(\  -  6)relgapx,  -3-2-1/2(5) 

Proof.  Let  H(e)  =  D{A  +  e  ■  (5/1/11(5/1112)1?.  Let  X,{()  be  the  t'-th  eigenvalue  of  H((), 
and  abbreviate  A,(0)  by  A,.  Let  relgap\,(e)  denote  the  relative  gap  of  the  i-ih  eigenvalue 
of  H((),  and  relgap\{a,b)  =  \a  ~  b\/{aby^^.  The  idea  is  that  if  €  is  small,  then  A,(e)  can 
only  change  by  a  small  relative  amount,  and  so  relgapx,{€)  can  only  change  by  a  small 
absolute  or  relative  amount.  Note  that  Aniin(/1)  can  decrease  by  as  much  as  ||(5,4||2.  Then 
by  Theorem  2.3  we  can  bound  relgap\,(()  below  by 

,  ,  ,  .     |A.(f)-A,(0|  ^     .     |A.-At|-(5(l-(5)-i(A. +  A,) 

relgapM     =     -- ^,,(,),^(,))./2  ^  ^I",?       (,.,,)i/2(i  ^  ,(i  _  ,)-.) 

>     (l-^)min(re/,ap.(A.,A,)-^.A_t0_) 

We  consider  two  cases,  relgap\(X,,  Xk)  >  2~^l^  and  Telgap\{Xi,Xif)  <  2~^i^.  The  first  case 
corresponds  to  A,  and  X^  differing  by  at  least  a  factor  of  2,  whence 

— ^ r;;  <  3  •  re/oapA(  A,,  At) 

(A,At)i/2  -  ^  ^^^         *' 

The  second  case  corresponds  to  A^  and  A^  differing  by  at  most  a  factor  of  2,  whence 

At    +   Afc         ^    ^       r.-l/2 


(A.A,)i/2 
Altogether  we  have 


<  3-2- 


/          3(5    \  /                    3-2-1/2  .^\ 
relgapxM  >  (1  -  ^)  (  1  -  ^Ts )  \^^^9apXi ^Z^ — I 

Now  integrate  the  bound  of  Theorem  2.7  from  e  =  0  to  e  =  ||5A||2  to  get  the  desired  result. 
I 

In  complete  analogy  to  [2],  we  may  also  prove 
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Proposition  2.10  Let  Aj   <  ■  •■  <  X„  be  the  eigenvalues  of  H  and  hi  <  ■■■  <  hn  be  its 

diagonal  entries  m  increasing  order.    Then 

h, 

In  other  words,   the  diagonal  entries  of  H  can  differ  from  the  eigenvalues  only  by  factors 
bounded  by  k(A). 

Proof.    See  the  proof  of  Proposition  2  in  [2]. 

Proposition  2.11   Let  H  -  DAD  with  eigenvalues  A,.  Let  d,  be  the  diagonal  entries  of  D. 
Let  u,  be  the  i-th  eigenvector  of  H  normalized  so  that  its  i-th  component  v,(i)  =  1.    Then 


1/2 


We  also  have 


\v,{j)\<v,U)  =  {K{A)f^'-min{lj-\        .  (^ 


\i\(j]\<{K(A]f^'-miR(j-,  ^) 

Uj      a, 


1/2 


In  other  words,  the  eigenvectors  are  scaled  analogously  to  the  diagonal  of  H . 

Proof.    See  the  proof  of  Proposition  6  in  [2]. 

Proposition  2.12   Let  H(e)  and  v^(e)  be  as  in  Theorem  2.1,  andVi(j)  be  as  in  Proposition 
2.11.    Then 


U'.(e)U)-t'.(0)(;)|< 


(2n- 2)1/2 


•^nun(A)-min(re/5apA,,2-i/2) 


€-U.(j)  +  0(€2) 


In  other  words,  each  component  of  each  eigenvector  is  perturbed  by  a  small  amount  relative 
to  its  upper  bound  of  Vi(j)  of  Proposition  2.11.  Thus  small  components  of  eigenvectors 
may  be  determined  with  as  much  relative  accuracy  as  large  components.  Note  that  relgapx, 
exceeds  2~'/^  only  when  A,  differs  from  its  nearest  neighbor  by  at  least  a  factor  of  2. 

Proof.    See  the  proof  of  Theorem  7  in  [2]. 

We  illustrate  these  results  with  two  examples.  First  we  consider  the  matrix  H  =  DAD 
of  the  introduction: 


H  = 


10"°  1029  10^^ 
1029  lo^o  10^ 
lO^^     10^        1 


,    A  = 


1  .1  .1 
.1  1  .1 
.1     .1     1 


and       D  =  diag(10'",10'M) 


To  six  correct  figures,  H's  eigenvalue  matrix  A  and  eigenvector  matrix  V  (normalized  to 
have  the  largest  entry  of  each  eigenvector  equal  to  1)  are 

A  =  diag  (1.00000-  10"°,  9.90000-  10*^,9.81818  ■  10"') 
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and 


V  = 


1.00000 
1.00000-  10-'^ 
1.00000-10-2^ 


-1.00000-  10-^' 

1.00000 
9.09091-  10-12 


-9.09091  •  10-" 

-9.09091  •  10-12 

1.00000 


One  may  compute  that  k{H)  ^  10"'°  and  /<(/!)=  1.33.  Thus,  according  to  Theorem  2.3. 
changing  each  entry  of  H  in  its  7th  decimal  place  or  beyond  would  not  change  A  in  the  figures 
shown.  The  refined  error  bounds  of  Proposition  2.5  are  essentially  the  same  in  this  case. 
One  can  further  verify  the  assertion  of  Proposition  2.10  that  the  ratios  of  the  eigenvalues  to 
the  diagonal  entries  of  H  are  bounded  between  .9  =  Amin(.4)  and  1.2  -  A^ax (-'!)•  One  may 
also  compute  that  the  relative  gaps  reigap\,  for  all  three  eigenvalues  are  appro.ximately  10^°. 
Thus,  according  to  Theorem  2.7.  7th  figure  changes  in  H  would  not  change  its  eigenvectors 
by  more  than  10-'^  in  norm.  In  fact,  the  eigenvectors  are  even  more  accurately  determined 
than  this.  Let  \'  =  {v,{j)}  be  the  matrix  of  upper  bounds  of  entries  of  V  as  define  in 
Proposition  2.11: 

1.5  1.5-  10-1°     1.5-  10-2°  " 

1.5-10-1°  1.5  1.5-10-1° 

1.5-10-2°     1.5.10-10     1.5-10-2° 


V'  =i 


Then  according  to  Proposition  2.12.  7th  figure  changes  in  H  cause  changes  in  at  most  the 
5th  digits  of  all  the  entries  of  V .  In  other  words,  for  this  examples  all  the  eigenvalues  and 
aJl  the  components  of  all  the  eigenvectors  are  determined  to  nearly  full  relative  precision 
by  the  data.  Later,  we  will  show  Jacob!  can  compute  them  with  this  accuracy.  In  contrast, 
QR  does  not  even  get  the  signs  of  the  two  small  eigenvalues  or  many  components  of  the 
eigenvectors  correct. 

The  second  example  serves  to  illustrate  the  difference  between  Theorem  2.3  and  the 
refined  bounds  of  Proposition  2.5.  Let  H  =  DAD  where  D  is  the  same  as  before  and 


A  = 


1        I-  fi     1  -M 

1-/2  1  1   —  M 

I-  fi        \-  fl  1 


where  fi  =  10-^.  The  eigenvalues  o{  H  are  10''°,  2-101''  and  1.5-10-^.  Now  k(A)  ss  10^,  so 
according  to  Theorem  2.3,  an  tj  relative  change  in  the  matrix  entries  will  cause  as  much  as 
a  lO^T]  relative  change  in  the  eigenvalues.  In  contrast,  the  refined  bounds  predict  a  relative 
change  of  rj  in  10''°  and  10^7/  in  the  two  smaller  eigenvalues.  Thus,  the  largest  eigenvalue  is 
just  as  insensitive  as  predicted  by  standard  norm  based  perturbations  theory. 

2.2      Optimality  of  the  Bounds  for  Symmetric  Positive  Definite  Matrices 

In  this  section  we  show  that  the  bounds  of  the  last  section  are  attainable.  In  particular, 
we  give  explicit  small  componentwise  relative  perturbations  which  attain  the  eigenvalue 
bounds;  this  will  imply  that  as  long  as  the  initial  matrix  entries  contain  small  relative  errors, 
the  eigenvalues  computed  by  Jacobi  (modulo  the  assumption  on  max,„  K(Am)/«(j4o);  see 
section  6)  or  bisection  are  optimally  accurate.  In  fact,  it  will  suffice  that  the  diagonal  entries 
of  matrix  have  small  relative  errors.  We  have  (necessarily)  slightly  weaker  results  for  the 
optimality  of  our  eigenvector  bounds. 
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We  begin  by  showing  that  the  assumption  \\SA\\2  <  Amin(  -•! )  of  the  last  section  is  essential 
to  having  relative  error  bounds  at  all.  If  this  bound  were  violated.  .4  +  6 A  (and  so  H  +  SH) 
could  become  indefinite,  implying  that  all  relative  accuracy  in  at  least  one  eigenvalue  is 
completely  lost.  In  contrast  to  standard  perturbation  theory,  however,  which  assumes 
a  bound  on  ||(*^//^|l2  instead  of  [[(5.4||2.  one  cannot  say  which  eigenvalue  will  lose  relative 
accuracy  first.  In  the  conventional  case,  as  ||(''/f||2  grows,  it  is  the  smallest  eigenvalues 
which  lose  accuracy  first,  the  larger  ones  remaining  accurate.  As  ||<?i.-l||2  grows,  however, 
any  eigenvalue  in  the  spectrum  (except  the  very  largest)  may  lose  its  relative  accuracy  first. 
The  following  example  illustrates  this: 


H  = 


10 


20 


1       .99 
.99      1 


10 


-20 


.-1  = 


1 


1 

.99 


.99 
1 


.  D  =  diag  (10'".  1.1,10 


i-10n 


Note  that  A„ijn(.-1)  =  .01.  .\s  \\SA\\2  approaches  .01,  the  eigenvalues  near  10^°.  1.99  and 
10~^°  retain  their  accuracy,  but  the  one  near  .01  can  lose  all  its  relative  accuracy. 

We  next  show  that  the  relative  error  bound  of  Theorem  2.1  can  be  nearly  attained  for 
at  least  one  eigenvalue  simply  by  making  appropriate  small  relative  perturbations  to  the 
diagonal  of  H . 

Proposition  2.13  Let  H  =  DAD  be  symmetric  positive  definite,  with  D  =  diag  (H^P) 
diagonal  and  A„  =  1.  Let  6A  =  rjl.  0  <  q  <  X^^i^),  arid  H  +  6H  =  D{A  +  6A)D.  Then 
for  some  i  we  have 


X,(H  +  SH] 


X,{H} 


>(!  + 


Amin(^"i) 


a/n 


Si   1  + 


nAniin(A) 


Proof.    We  have 


YlX.iH)  =  det(DAD)  =  det(D^)det(A)  =  det(Z}2)  JJ  A,(A) 


and 


Y[X,{n  +  6H)  =  det{D(A  +  r,I)D)  =  det(Z)2)  det(A  +  r?/)  =  dei{D'^)Y[[X,[A)  +  /?) 


Therefore 


y.X,{H  +  bH)  ^.p^A.(A)  +  r? 


K{H) 


K{A) 


Aminl  -4 ) 


implying  that  at  least  one  term  A,(/r  +  6H)/Xf{n)  must  exceed  (1  +  T?/Anun(^))^''"-  This 
last  expression  is  approximately  1  +  r?/(nAnun(>l))  when  rj  C  Ainjn(^)-       I 

The  example  at  the  beginning  of  this  section  showed  that  the  error  bound  of  Theorem 
2.3  and  the  last  proposition  may  only  be  attained  for  one  eigenvalue.  Proposition  2.5  of  the 
last  subsection  showed  that  for  asymptotically  small  ||<5/1||2,  the  maximum  perturbation  in 
each  eigenvalue  may  be  attained  only  with  small  diagonal  perturbations  of  A. 
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After  we  show  that  the  rounding  errors  introduced  by  Jacobi  are  of  the  form  ||(''.4||2  = 
0(e)  in  the  next  section.  Propositions  2.13  and  2.5  will  show  that  Jacobi  (modulo  the 
assumption  on  max^  «(-•!„)/«(  .4o))  computes  all  the  eigenvalues  will  optimal  accuracy, 
provided  that  only  the  diagonal  entries  of  H  have  small  relative  errors.  The  same  optimality 
property  is  true  of  bisection. 

Now  we  consider  eigenvectors.  Here  our  results  are  necessarily  weaker,  as  the  following 
example  shows.  Suppose  H  is  diagonal  with  distinct  eigenvalues.  Then  small  relative 
perturbations  to  the  matrix  entries  leave  H  diagonal  and  its  eigenvalue  matrix  (the  identity 
matrix)  unchanged.  Therefore,  the  only  way  we  can  hope  to  attain  the  bounds  of  Theorem 
2.7  is  to  use  perturbations  6 A  which  are  possibly  dense,  even  if  H  is  not.  Furthermore,  a 
block  diagonal  example  like  the  first  one  in  this  section  shows  that  the  attainable  eigenvector 
perturbations  will  not  necessarily  grow  with  k(.-I).  Thus,  the  best  we  can  prove  is 

Proposition  2.14  Let  H  =  DAD.  A,,  i',,  SH,  6A  and  t]  be  as  in  Proposition  2.5.  Let  i-[ 
be  the  unit  eigenvectors  of  H  -\-  811 .  Then  one  can  choose  SA,  ||<5.4||2  =  77  <C  Ajnin(^4),  50 
that 

lit'.  -  v:h  >  T -rfl +  O(v') 

Proof.  Consider  expression  (2.8)  for  v,  —  v[  (there  8A  is  written  (E).  By  using  a 
Householder  transformation,  one  can  prove  there  exists  a  symmetric  6  A  such  that  y^SAy,  = 
||2/fc||2||y.|l2ll<5-4||2  for  arbitrary  yk  and  j/,.  Since  A^  =  yl Ayk  <  ||3/*ipAnu»x(^).  we  can  find 
6A  to  make  yJiAy,  >  (A,Afc)^''^||(!)/l||2/Aniax(^)-  Choosing  k  so  that  Afc  is  closest  to  A, 
completes  the  proof.       I 

2.3      General  Matrices 

The  results  on  singular  values  and  singular  vectors  are  analogous  to  the  results  for  eigen- 
values and  eigenvectors  in  the  first  subsection.  Just  as  we  derived  perturbation  bounds  for 
eigenvalues  from  a  more  general  result  for  generalized  eigenvalues  of  pencils,  we  will  start 
with  a  perturbation  bound  for  generalized  singular  values  and  then  specialize  to  standard 
singular  values. 

Let  G\  and  G2  be  matrices  with  the  same  number  of  columns.  G2  of  full  column  rank, 
and  otherwise  both  arbitrary.  We  define  the  i-th  generalized  singular  value  <7,(Gi,  G2)  of  the 
pair  (Gi,  G2)  as  the  square  root  of  the  i-th  eigenvalue  of  the  definite  pencil  G^G\  -  XG2G2 
[9].  If  we  let  G2  be  the  identity,  ct,{Gi,  G2)  is  the  same  as  the  standard  singular  value  (7,(Gi ) 
ofGi. 

Lemma  2.15  Let  Gi  and  G2  be  matrices  with  the  same  number  of  columns,  G2  of  full 
column  rank,  and  otherwise  both  arbitrary.  Let  6Gj  be  a  perturbation  of  Gj  such  that 

WiGjxh  <  ijWGjxh 

for  all  X  and  some  rj^  <  1.  Let  Oi  be  the  i-th  generalized  singular  value  of  (G\,G2)  and  a[ 
be  the  i-th  generalized  singular  value  of  [Gi  +  6Gi,G2  +  <5G2).   Then  either  Oi  =  a,'  =  0  or 

1  -  ^1  <  ^  <-  1  +  ^1 
1  -I-  772  ~  a,  ~  1  -  7?2 
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Proof.    From  the  Courant-Fischer  minimcLX  theorem  [12]  we  have 

'\Gil\\2 


a,  —  mm  max 

S'   xeS'  II^2-e||2 

where  the  minimum  is  over  all  i-dimensional  subspaces  S'.  The  rest  of  the  proof  is  analogous 
to  that  of  Lemma  2.1.      I 

Lemma  2.16  Let  G\  and  G2  be  as  in  Lemma  2.15.  Let  G_,  =  Bj^^  where  Aj  has  full  rank 
and  is  otherwise  arbitrary.  Let  6Gj  =  SBjAj  be  a  perturbation  of  Gj  such  that  \\bBjX\\2  < 
T]j\\Bjx\\2  for  all  x  and  some  r/j  <  1.  Let  a,  and  a\  be  the  i-th  generalized  singular  values  of 
{G\,G2)  and  (Gi  +  6Gi,G2  +  SG2).  respectively.   Then  either  a,  =  a[  =  0  or 

1  -  ^1  <-  ^  <  1  +  ^1 
1  +  T/2  ~  cr,  -  1  -  r?2 

The  proof  is  analogous  the  proof  of  Lemma  2.2. 

Theorem  2.17  Let  G  -  BD  be  a  general  full  rank  matrix,  and  D  chosen  diagonal  so  that 
the  columns  of  B  have  unit  two-norm  (i.e.  D„  equals  the  two-norm  of  the  i-th  column  of 
G).  Let  6G  =  SBD  be  a  perturbation  of  G  such  that  \\SB\\2  =  77  <  (T„un{B).  Let  ct,  and  a[ 
be  the  i-th  singular  values  of  G  and  G  +  6G ,  respectively.   Then 

^■^^^<-^<^(B)-v  (2.18) 

where  n{B)  =  crmax(-S)/«Tmin(-5)  <  Ti^^^/<^mlniB),  and  n  is  the  number  of  columns  of  G.  In 
particular,  if  \6G,j/iG,j\  <  T]/n,  then  ||<5J9||2  <  t?  and  the  bound  (2.18)  applies. 

The  proof  is  analogous  to  that  of  Theorem  2.3. 

Just  as  the  bounds  of  Theorem  2.3  were  not  attainable  by  all  eigenvalues,  neither  are 
the  bounds  of  Theorem  2.17  attainable  for  all  singular  values.  Antilogous  to  Proposition 
2.5,  we  may  derive  tighter  bounds  for  individual  singular  values. 

Proposition  2.19  Let  G  =  BD  be  as  in  Theorem  2.17,  with  singular  values  a,,  right  unit 
singular  vectors  Vi  and  left  unit  singular  vectors  Ui.  Let  G  ->r6G  =  {B  -\-  SB)D  have  singular 
values  a[,  where  ||^B||2  =  ^  <C  o'min(-B)-    Then  the  bound 

— '-  <  — —  +  0{r)^)  (2.20) 


is  attainable  by  the  perturbation  6B  =  rjUii^Dvi)^ /\\Dv,\\2. 

Proof.  Bound  (2.20)  is  derived  from  the  standard  first  order  perturbation  theory  which 
says  a,{G  +  6G)  =  a,(G)  +  uJSGv,  +  0(||<5G||^),  and  substituting  \uj6Gv,\  =  \uj6BDv,\  < 
\\Dv,\\2\\bB\\2.  This  last  inequality  is  dearly  attained  by  the  choice  of  SB  in  the  statement 
of  the  proposition.      I 
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Now  we  consider  the  singular  vectors.  For  simplicity  w'e  assume  G  is  square.  We  use  the 

'  V      V    ' 
U     -U 


fact  that  if  G  =  L"£\'^  is  the  singular  value  decomposition  of  G.  then  2   ^Z- 


is  the  eigenvector  matrix  of  the  svmmetric  matrix 


[9].    Therefore  we  can  use 


0     G^ 
G      0 

perturbation  theory  for  eigenvectors  of  symmetric  matrices  to  do  perturbation  theory  for 
singular  vectors  of  general  matrices. 

We  also  need  to  define  the  gaps  for  the  singular  vector  problem.  The  absolute  gap  for 
singular  values   is 

,  .    kt-(7fc| 

absgapo^  =  min  — -r-— — 

i.e.  essentially  the  same  as  the  absolute  gap  for  eigenvalues.   However  the  relative  gap  for 
singular  values 

,  ■   k.  -  <^k\ 

relgapff,  =  mm 

;c#.    a,  +  cTk 

is  somewhat  different  from  the  relative  gap  for  eigenvalues. 

The  standard  perturbation  theorem  for  singular  vectors  is  essentially  the  same  as  for 
eigenvectors.  Let  G  have  right  (or  left)  unit  singular  vectors  i',,  and  let  G  +  6G  have  right 
(or  left)  unit  singular  vectors  v[.  Let  rj  =  pG||2/||G||2-  Then 


k<  -  ^'ih  < 


absgapa^ 


+  0('?' 


We  improve  this  as  follows: 


Theorem  2.21  Let  G  =  BD  be  as  in  Theorem  2.17.  Define  G{e)  =  {B  +  eE)D  where  E 
is  any  matrix  with  unit  two-norm.  Let  (7,(e)  be  the  i-th  singular  value  of  G(e),  and  assume 
<7,(0)  is  simple  so  that  the  corresponding  right  unit  singular  vector  v,{e)  and  left  unit  singular 
vector  u,{()  are  well  defined  for  sufficiently  small  e.   Then 


max(||i',(e)-j;.(0)||2,  ||u.(e)-u.(0)||2)  < 


(n-.5)i/2f 


■+0(€2)< 


(n-  .5)i/2^(5)e 


(^miniB)  ■  relgap^,  '  ~''  '  relgap^, 

Proof.    Let  Vk{0)  be  abbreviated  by  v^  and  Ufc(O)  by  u^.  Define 


+  0(€^ 


G  = 


0     G^ 
G      0 


,  E  = 


0     E^ 
E      0 


,  D  = 


D    0 
0     / 


,  B  = 


0     B^ 
B      0 


,-f  = 


J_ 

7! 


±u, 


so  that  G  =  DBD.  G  has  eigenvalues  ±a,  corresponding  to  eigenvectors  i*.  Let  xf  (e)  be 
the  eigenvectors  of 


0  D(B^  +  €E^) 

iB  +  (E)D  0 


=  D{B  +  eE)b 


Then  from  [9]  we  have 


±i/±. 


±<^,  T  CTk 


xf  +  0{e') 


(2.22) 
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Now  we  compute 

\\'e  also  have 

so 

Taking  norms  in  (2.22)  yields  the  result.      I 

Corollary  2.23   Itt  G  -  BD  he  as  in  Theorem  2.17.  Suppose  6  =  \\SB\\2/(7nMn{B)  satisfies 

S 

<  relgap^, 


1  -6 

Let  u,  and  u,  be  the  unit  right  and  left  singular  vectors  of  G,  respectively,  and  let  v'  and  u' 
be  the  unit  right  and  left  singular  vectors  of  G'  =  {B  +  SB)D.  respectively.   Then 

(n  -   5)1/2(5 
majc(||t;,  -  t;,'||2.  ||u,  -  u[\\2)  < 


{l-6){(l-i)relgap^,-6) 


Proof.    The  proof  is  analogous  to  the  proof  of  Corollary  2.9.      I 
There  are  analogs  to  Propositions  2.10  through  2.12  of  the  last  section,  gotten  by  con- 
sidering H  =  G'^G: 

Proposition  2.24  Let  G  =  BD  be  as  in  Theorem  2.17.  Let  Oi  <■■■  <  a^  be  the  singular 
values  of  G  and  d\  <  ■  ■  ■  <  d^  the  diagonal  entries  of  D  in  increasing  order.    Then 

a, 

<^nua{B)  <    -J-  <  Crmax(5) 

d, 

Proposition  2.25  Let  G  =  BD  be  as  in  Theorem  2.17  with  singular  values  ctj  <  •  •  •  <  cr„. 
Let  v,  be  the  i-th  right  singular  vector  of  G,  normalized  so  that  its  i-th  component  r,(j)  =  1. 
Then 


|v,(j)|<5.0)  =  («(5))3.min(^,^) 


We  also  have 


\3       •   /  '^t    dj 


\v,iJ)\<{K{B)f-nnn(-^,-^) 

dj    a, 


Proposition  2.26  Let  G{€)  and  Viie)  be  as  in  Theorem  2.21,  and  v,{j)  as  in  Proposition 
2.25.    Then 

\vM){J)  -  vm{j)\  <     2  ^)"  "  ^^!^'        •  e  •  Hi)  +  0(6^) 
'^imniB)  ■  relgap^. 

There  are  analogs  to  ail  the  results  in  this  section  for  matrices  G  =  DB  scaled  from  the 
left  instead  of  the  right.  Thus,  one  can  choose  to  scale  either  the  rows  or  the  columns  of  G 
to  have  unit  two-norms,  whichever  one  minimizes  the  condition  number. 
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2.4      Optimality  of  the  Bounds  for  General  Matrices 

The  results  in  this  section  are  andogous  to  but  necessarily  weaker  than  the  results  of 
subsection  2.2.  In  particular,  it  is  no  longer  the  case  that  the  perturbation  bounds  for  the 
singular  \-alues  can  be  attained  by  small  relative  perturbations  in  the  matrix  entries. 

First  consider  the  restriction  \\SB\\2  <  (^nuniB).  Just  as  in  the  symmetric  positive 
definite  case,  this  is  necessary  so  that  B  +  SB  remains  nonsingular.  When  B  +  6B  becomes 
singular,  at  least  one  singular  value  will  necessarily  lose  all  relative  accuracy.  The  same 
l<ind  of  block  diagonal  example  as  in  subsection  2.2  also  shows  that  only  one  singular  value 
may  have  its  sensitivity  depend  on  k(B),  and  it  might  be  anywhere  in  the  spectrum  (except 
the  very  largest  singular  value). 

In  order  to  prove  an  analogue  of  Proposition  2.13,  we  must  permit  perturbations  SB  o(  B 
which  are  small  in  norm  but  may  make  large  relative  changes  in  tiny  entries  of  5  (a  similar 
perturbation  was  needed  to  prove  that  the  bound  in  Proposition  2.19  was  attainable): 

Proposition  2.27  Let  G  =  BD  with  D  diagonal  and  the  columns  of  B  having  unit  norm. 
Then  there  exists  a  SB  with  \\SB\\2  =  n  <  <7min(5)  such  that  for  G  +  SG  =  (B  +  SB)D  we 
have  for  at  least  one  i 

aAG  +  SG)^  V        ,;„^^^         V 


If  we  restrict  SB  so  that  |^5,j/5,j|  <  tj.  then  such  a  perturbation  SB  may  not  exist. 

Proof.  The  proof  is  very  similar  to  that  of  Proposition  2.13.  Let  A'  be  a  rank  one 
matrix  of  minimal  2-norm  such  that  B  +  A'  is  singular,  and  let  SB  =  —tjX.  Then  as  in 
Proposition  2.13  we  discover  that 

yra,{G  +  SG)  ^         _J]__ 

V         '^'iG)  (^rrun(B) 

and  so  at  least  one  term  (t,{G  +  SG)/a,(G)  exceeds  (1  +  rj/'^miniB))^^^.  To  see  that  small 
componentwise  relative  perturbations  are  not  sufficient,  consider  the  matrix 

1      1 


G  =  B  = 


—e     € 


with  e  C  1.  The  condition  number  of  B  is  approximately  l/e,  and  relative  perturbations 
of  size  T]  in  its  entries  cannot  change  its  singular  values  by  more  than  a  factor  of  about 

As  in  Proposition  2.14,  our  lower  bound  on  the  attainable  perturbations  in  the  singular 
vectors  requires  a  dense  SB  and  does  not  grow  with  k(B). 

Proposition  2.28  Let  G  =  BD,  a,,  u,,  r,,  SG  =  SBD  and  tj  be  as  in  Proposition  2.19. 
Let  u[  and  v[  be  the  unit  left  and  right  singular  vectors  of  G  +  SG,  respectively.  Then  one 
can  choose  SB,  ||^5||2  =  t?  <  (yrtaa{B),  so  that 

max(||u.  -  u:||2,  lit;.  -  v[h)  >  ,L,,,,,„     +  0{v') 

^^'^CTiii^iBjrelgapg, 

Proof.  Consider  expression  (2.22)  for  the  perturbation  (ti,  -  u[,Vi  -  v[)  (there  SB  is 
written  eE).  Choose  k  so  jcr,  -  a^l  is  minimized.  Let  SB  =  T]Ui{Dvk)^ /\\Dvk\\2  if  c^^fc  >  cr^ 
and  SB  -  77Ufc(L't;,)^/||£)u,||2  otherwise.  Now  compute.      I 
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3      Two-sided  Jacobi 

In  this  section  we  prove  that  two-sided  Jacobi  in  floating  point  arithmetic  applied  to  a 
positive  definite  symmetric  matrix  computes  the  eigenvalues  and  eigenvectors  with  the  error 
bounds  of  section  2. 

In  this  introduction  we  present  the  algorithm  and  our  model  of  floating  point  arithmetic. 
In  subsection  3.1,  we  derive  error  bounds  for  the  computed  eigenvalues.  In  subsection  3.2, 
we  derive  error  bounds  for  the  computed  eigenvectors. 

Let  Hq  -  DqAoDq  be  the  initial  matrix,  and  Hm  =  DmAmD-m  where  Hm  is  obtained 
from  Hm-\  by  applying  a  single  Jacobi  rotation.  Here  Dm  is  diagonal  and  Am  has  unit 
diagonal  as  before.  All  the  error  bounds  in  this  section  contain  the  factor  max^^j  k(  .-!„ ), 
whereas  the  perturbation  bounds  of  section  2  are  proportional  to  k{Aq).  Therefore,  our 
claim  that  Jacobi  solves  the  eigenproblem  as  accurately  as  predicted  in  section  2  depends 
on  the  ratio  max„,  K(.4m  )/«(^o)  being  modest  in  size.  Note  that  convergence  of  Hm  to 
diagonal  form  is  equivalent  to  the  convergence  of  Am  to  the  identity,  or  K(Am]  to  1.  Thus 
we  expect  K(.4m)  to  be  less  than  k(.4o)  eventually. 

We  have  overwhelming  numerical  evidence  that  maXm  n{Am)l >^[Aq)  is  modest  in  size:  in 
section  7.  the  largest  value  this  ratio  attained  was  1.82  in  over  900000  numerical  experiments. 
Our  theoretical  understanding  of  why  this  ratio  is  so  small  is  somewhat  weaker;  we  present 
our  theoretical  bounds  on  this  ratio  in  section  6. 

The  essential  difference  between  our  algorithm  and  standard  two-sided  Jacobi  is  the 
stopping  criterion;  according  to  Theorem  2.3,  we  must  set  H^j  to  zero  only  li Hij/iHuHjj)^^'^ 
is  small,  not  just  if  //,j/mcLXfc/  \Ht;i\  is  small.  This  stopping  criterion  has  been  suggested 
before  [16,5,3],  but  without  our  explanation  of  its  benefits.  Otherwise,  our  algorithm  is  a 
simplification  of  the  standard  one  introduced  by  Rutishauser  [12].  We  have  chosen  a  simple 
version  of  the  algorithm,  omitting  enhancements  like  delayed  updates  of  the  diagonals  and 
fast  rotations,  to  make  the  error  analysis  clearer  (an  error  analysis  of  these  enhancements 
is  future  work). 

Algorithm  3.1  Two-sided  Jacobi  for  the  symmetric  positive  definite  eigenproblem.  tol  is  a 
user  defined  stopping  criterion.  The  matrix  V  whose  columns  are  the  computed  eigenvectors 
initially  contains  the  identity. 


repeat 

for  all  pairs  i  <  j 


/*  compute  the  Jacobi  rotation  which  diagonalizes 

C  =  (6-a)/(2c) 


Hit     Hij 


a     c 

c     b 


r 


<  =  5i^n(c)/(ici  +  yr+c^) 
cs  =  i/yrn^ 

sn  =  cs  *  t 

/*  update  the  2  by  2  submatrix  */ 

Hi,  —  a  -  c  *  t 

Hjj  =  b  +  c*t 

H,j  =  H„  =  0 

/*  update  the  rest  of  rows  and  columns  i  and  j  */ 
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for  k  —  I  to  n  except  i  and  j 
imp  =  H^k 

Htk  =  cs  *  tmp  -  sn  *  Hjk 
Hjk  =  sn  *  tmp  +  cs  *  Hjk 

Hk:   =   H,k 

Hkj  =  Hjk 
endfor 

/*  update  the  eigenvector  matrix  V  */ 
for  k  =  I  to  n 

tmp  =  Vk, 

Vk,  =  cs  *  tmp  -  sn  *  Vkj 

Vfcj  =  5n  ♦  tmp  +  cs  *  Vkj 
endfor 
endfor 
until  convergence  (all  \H,j\/{H„Hjj)^^^  <  tol  ) 

Our  model  of  arithmetic  is  a  variation  on  the  standard  one:   the  floating  point  result 
//(■)  of  the  operation  (•)  is  given  by 


(3.i; 


where  |£,|  <  s,  and  £  <C  1  is  the  machine  precision.  This  is  somewhat  more  general  than  the 
usual  model  which  uses  fl{a±h)  =  {a±h){l-\-£i)  and  includes  machines  like  the  Cray  which 
do  not  have  a  guard  digit.  This  does  not  greatly  complicate  the  error  analysis,  but  it  is 
possible  that  the  computed  rotation  angle  may  be  off  by  a  factor  of  2,  whereas  with  a  guard 
digit  the  rotation  angle  is  always  highly  accurate.  This  may  adversely  affect  convergence, 
but  as  we  will  see  it  does  not  affect  the  one-step  error  analysis. 

Numerically  subscripted  es  will  denote  independent  quantities  bounded  in  magnitude 
by  £.    As  usual,  we  will  make  approximations  like  (1  +  i£i){l  +  J£2)  =  1  +  («  +  j)sz  and 

(l  +  t£l)/(l  +  J£2)=   l  +  (l  +  ;)£3. 

3.1      Error  Bounds  for  Eigenvalues  Computed  by  Two-sided  Jacobi 

The  next  theorem  and  its  corolljiry  justify  our  accuracy  claims  for  eigenvalues  computed  by 
two-sided  Jacobi. 

Theorem  3.2  Let  Hm  b^  the  sequence  of  matrices  generated  by  Algorithm  3.1  in  finite 
precision  arithmetic  with  precision  e;  that  is  Hm+\  is  obtained  from  Hm  by  applying  a 
single  Jacobi  rotation.   Then  the  following  diagram  commutes. 

U     floating,     rr 
■"m  Jacob;  V-"m -1-1 

^6H. 


fl{a±b) 

=     a(l +£i)±6(l-f  £2) 

fliaxb) 

=     (ax6)(l  +  £3) 

fl{a/b) 

=     (a/6)(l  +  £4) 

//(v^) 

=       ^Ail+e5) 
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The  top  arrow  indicates  that  Hm  +  i  is  obtained  from  Hm  by  applying  one  Jacohi  rotation 
in  floating  point  arithmetic.  The  diagonal  arrow  indicates  that  Bm+i  '■?  obtained  from  H'^  by 
applying  one  Jacobi  rotation  in  exact  arithmetic;  thus  Hm  +  i  a'^'^  ^m  ^^^  exactly  similar. 
The  vertical  arrow  indicates  that  H'^  =  Hm  +  ^Hm-  ^Hm  is  bounded  as  follows.  Write 
SHm  =  Pm^AmDm-    Then 


\6Amh  <  (182(2n- 4)^/2 +  73)£ 


(3.3) 


In  other  words,  one  step  of  Jacobi  satisfies  the  assumptions  needed  for  the  error  bounds  of 
section  2. 

Corollary  3.4  Assume  Algorithm  3.1  converges,  and  that  H\f  is  the  final  matrix  whose 
diagonal  entries  we  take  as  the  eigenvalues.  Write  Hm  —  DmAmDm  with  Dm  diagonal  and 
Am  '^'ith  ones  on  the  diagonal  for  0  <  m  <  M .  Let  \j  be  the  j-th  eigenvalue  of  Hq  and  A' 
be  the  j-th  diagonal  entry  of  Hm.    Then  to  first  order  in  £  the  following  error  bound  holds: 


lA,  -  A' 


A, 


•"<  (5 -M -(182(271-4)^/^ +  73) +  n-<o/)-    max    K(Ar 

0<m<A/ 


(3.5) 


Remark.  In  numericaJ  experiments  presented  in  section  7,  there  was  no  evidence  that  the 
error  bound  in  (3.5)  grew_with  increasing  n  or  M. 

Proof  of  Corollary  3.4.  Bound  (3.5)  follows  by  substituting  the  bound  (3.3)  and 
the  stopping  criterion  into  Theorem  2.3.       I 
Remark.  A  similar  bound  can  be  obtained  based  on  the  error  bound  in  Proposition  2.5. 

Proof  of  Theorem  3.2.      The  proof  of  the  commuting  diagram  is  a  tedious  compu- 
tation. Write  the  2  by  2  submatrix  of  Hmm  being  reduced  as 


a    c 

c     b 


d^       zdfdj 
zd,dj       d^ 


where  we  assume  without  loss  of  generality  that  a  >  b  and  c  >  0.  By  positive  definiteness 
0  <  z  <  z  =  (K{Am)  -  l)/{K{Am)  +  1)  <  1.  Let  a'  and  b'  be  the  new  values  of  H,,  and 
Hjj  computed  by  the  aJgorithm,  respectively.  Let  x  =  d^/df  <  1.  We  consider  two  cases, 
X  <  X  =  (v/5  -  l)/2  a:  .62,  and  x  >  x. 

First  consider'i  <  x.  Systematic  application  of  formulas  (3.1)  shows  that 

C     =     A(6-a)/(2*c)) 

=     (1  +  f4)(((l  +  £i)6  -  (1  +  £2)a)/{{l  +  e3)2c)) 
(l  +  £4)(l  +  £2)  f'b-a' 


l  +  £3 


2c 


where  6  =  (1  +  £i)b/{l  +  £2)  =  (1  +  £b)b,  \£b\  <  2£.  Thus  C  =  (1  +  £c)(^  -  a)/(2c)  where 

Let  tic)  denote  the  true  value  of  (  (i.e.  without  rounding  error)  as  a  function  of  a,  6  and 
c.  Using  (3.1)  again  one  can  show  i  =  (1  +  St)t{c)  where  \et\  <  7£. 
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Next 


b'    =     fl(b+ct)^{\  +  ei)b  +  (l  +  se)(l  +  £7)ct 

(oH -, :7-, ; c*(c) 


fine  the  exact  Jacobi  rotation  J™   = 


1  +  fl  (l  +  52)(l  +  >^5; 

=     (l  +  Sb'){b+{l  +  s,t(c))ct{c))  (3.6) 

where  l^ctjc)!  <  125  and  \sb'\  <  35.   Since  |^(c)|  is  an  increasing  function  of  c,  we  can  write 
(1  +£ct{c))cnc)  =  (l+Sc)c-t{(l  +  £jc),for  some  5,  where  \Sc\  <  \£ct(c)\  <  125. 

Now  we  can  define  c  =  ( 1  +  5c)c.  and  (,  i.  cs  and  s'n  as  the  true  values  of  the  untilded 
quantities  computed  without  rounding  error  starting  from  a,  b  and  c.    cs  and  s'n  will  de- 

CS         S  Tl 

which  transforms  H'    to  i/m+i  in  the 
—  sn     cs 

commutative  diagram  in  the  statement  of  the  theorem:  J^H'^Jm  —  Hm+i- 

Now  we  begin  constructing  SHm-  ^Hm  will  be  nonzero  only  in  rows  and  columns  i  and 
j.  First  we  compute  its  entries  outside  the  2  by  2  {i,j)  submatrix.  Using  (3.1)  one  can 
show  cs  =  (1  +  5cj)c~5  and  sn  =  {I  +  Ssn}sn  where  \Scs\  <  225  and  \£,n\  <  305.  Now  let  H',!^ 
and  H'  1^  denote  the  updated  quantities  computed  by  the  algorithm..  Then 

^'ik     =  fl{cs*  H,k  -  sn*  Hjk) 

=  ( 1  +  5io)(  1  +  58 )csH,k  -  ( 1  +  59)(  1  +  5n  )snn,k 

=  (1  +  5io)(l  +  £s){'i-  +  Sc,)csH,k  -  (1  +  59)(1  +  5n)(l  +  £,n)snHjk 

=  csH.k-  snHjk  +  i{H[k)  (3.7) 

Similarly 

H'jk     =     fl{sn*n,k  +  cs*H,k) 

=      (1  +  5i4)(l  +  5i2)(l  +  £sn)snH,k  +  (1  +  £l3)(l  +  £l5)(l  +  £c,)csHjk 

=    snn,k  +  csHjk  +  €{H'^,,)  (3.8) 

Now  X  =  dj/d,  implies 

-b-a^d]-d^^x^-l 

2c     ~   2zd,dj    ~     2zx 
where  i  =  2(1  +  £c)-  Then  x  <  x  implies 

1  zx 


t   = 


^+(i-(w)0"^"'"^ 


Also  |5n|  <  |f|,  so  this  last  expression  is  an  upper  bound  on  \s'n\  as  well.  Substituting  this 
bound  on  s'n,  cs  <  1,  1.^,^1  <  d^dkZ  and  \H-jk\  <  djdkz  into  (3.7)  and  (3.8)  yields 

\({H',k)\     <     56£d,dkz 
\e{H'^k)\     <     56£d,dk-z/(l-x^) 
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Thus 


Jin 
Jl 


Jl- 


//,fc  1 

■    f(//, 

k)  ' 

\ii>k  ^ 

T 

_  i(Hjk)  _ 

1 
\ 

+  Jm   ■ 

/ 
\ 

-I- 

SH 

,k 
J''  . 

J 

where 


\8Hjk\     <     n2£djdkz/{l-  x'^] 


Now  we  construct  the  2  by  2  submatrix  A  o{  SHm  at  the  intersection  of  rows  and  columns 
i  and  j.  We  will  construct  it  of  three  components  A  =  Ai  +  Aj  +  A3. 

Consider  the  formula  a'  =  fl(a  —  c  *  t)  for  the  i,  i  entry  of  Hm+\-  Applying  (3.1) 
systematically,  we  see 

a'     =     (l  +  eis)a-{l  +  eiT){l  +  £i6)ct 

=     (1  ^  £is)a  -  (1  +  £i7)(l  +  f i6)(l  +  £t]ct(c) 

,,   ^  ,  (l+gl7)(l  +  £l6)(l+g,)ct(c) 

=     (l  +  £i8)a — 

=     (l  +  £i8)a-(l  +  4(,))ct(c) 

where  kct(c)l  -  ^l^-  Since  a  >  0  and  ct(c)  <  0,  we  get 

figa  -  £'  ,  vCt(c)\ 
a'={l  + '^^ (a  -  ct(c))  =  (1  +  e,.){a  -  ct(c)) 


a  —  ct{c) 


where  |fa'l  <  2U.* 
Now  let 


Ai  = 


From  earlier  discussion  we  see 


a    c 

c     b 


0        £cC 

£cC      £bb 


+A^]Jm= 


0  C-  C 

c - c    6-6 


a  -  ct{c)  0 

0  b  +  ct{c) 


Next  let 


A2    =  £a' 


a     c 

c     b 


+  Ai 
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Thus 


Ji 


a     c 
c     b 


+  Ai   +  A2       Jm    =  il  +  Sa' 


a  -  ct{c)  0 

0  b  +  ct(c) 


a'  0 

0       67((1  +  £,0(1  +  ^6) 


Finally,  let 

A3  =   Jn 


0  0 

0      6'(1-1/((1  +£,0(1  +  ^6)) 


J^  = 

•J  m     — 


where  |56"|  <  k^'l  +  k6l  <  24e.  Then 


a     c 

c     h 


+  Ai   +  A2  +  A3   Um   = 


a'     0 
0     b' 


as  desired.  This  completes  the  construction  of  SHm-  We  may  bound 

112(2n- 4)1/^5 


11^.4 


m||2   < 


1-i^ 


+  73    £ 


(3.9) 


Now  we  consider  the  second  case,  when  x  >  x.  The  only  thing  that  changes  in  the 
previous  analysis  is  our  anjilysis  of  6H,k  and  6Hjk,  since  s'n  is  no  longer  small.  Instead  we 
substitute  the  bounds  \s'n\  <  1,  \cs\  <  1,  \H,it\  <  d,dkZ  <  djdkz/x  and  \Hjk\  <  djd^z  into 
(3.7)  and  (3.8)  to  get 

|€(//:^)|       <       b6£d,dkZ 

\i{H'jk)\     <     b^ed^dkz 
whence 


and 


\\6A 


\SH,k\ 
\^H,k\ 

m||2<    I 


<  U2ed,dkz 

<  I12£djdkz/x 


'112(271-4)1/22 


+  73    £ 


(3.10) 


Fincilly,  we  note  our  choice  of  x  makes  the  upper  bounds  in  (3.9)  and  (3.10)  both  equal, 
with  1/(1  -  i^)  =  1/i  <  1.62,  proving  the  theorem.       I 

Remark.  The  quantity  182(2n-4)'/2  in  the  theorem  may  be  multiplied  by  majCm,,^j  l^m.^l  < 
1.  Thus  if  the  Am  are  strongly  diagonally  dominant,  the  part  of  the  error  term  which  de- 
pends on  n  is  suppressed. 

Commutative  diagrams  like  the  one  in  the  theorem,  where  performing  one  step  of  the 
algorithm  in  floating  point  arithmetic  is  equivalent  to  making  small  relative  errors  in  the 
matrix  and  then  performing  the  algorithm  exaictly,  occur  elsewhere  in  numerical  analysis. 
For  example,  such  a  diagram  describes  an  entire  sweep  of  the  zero-shift  bidiagonal  QR 
algorithm  [7],  and  is  the  key  to  the  high  accuracy  achieved  by  that  algorithm.  Also,  if  one 
modifies  one  line  in  the  traditional  zero-shift  symmetric  tridiagonal  QR  algorithm  so  it  is 
performed  in  higher  precision,  such  a  diagram  describes  an  entire  sweep  of  that  algorithm  as 
well  [11].  Unfortunately,  this  ability  to  have  a  high  accuracy  sweep  does  not  often  translate 
into  overall  accuracy,  because  in  tridiagonal  QR  k(Ato)  frequently  grows  greatly  for  that 
algorithm  before  converging  to  1. 
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3.2      Error  Bounds  for  Eigenvectors  Computed  by  Two-sided  Jacobi 

The  next  two  theorems  justify  our  accuracy  claims  for  eigenvectors  computed  by  two-sided 
Jacobi. 

Theorem  3.11  Let  V  =  [I'l iv]  be  the  matrix  of  unit  eigenvectors  computed  by  Algo- 
rithm 3.1  in  finite  precision  arithmetic  with  precision  s.  Let  Vj  =  [vji, ....  vj^]  be  the  true 
eigenvector  matrix.  Let  R  =  rmLKm  K{Am)  be  the  largest  K(Am)  of  any  iterate.  Then  the 
error  in  the  computed  eigenvectors  is  bounded  in  norm  by 

n  II     ^  ("  -  l)^/'(n  ■?o/  +  .U  .(182(271-4)^/    +73)6)k   ,    ^^,, 

K'.  -  i'T<   2  <  ; +  46A/5  (3.12) 

relgaps. 

Remark  In  numerical  experiments  presented  in  section  7,  there  was  no  evidence  that  the 
error  bound  in  (3.12)  grew  with  increasing  n  or  M. 

Proof.   Let  Hq H\r  be  the  sequence  of  matrices  generated  by  the  Jacobi  algorithm. 

where  H\{  satisfies  the  stopping  criterion.  Let  Jm  be  the  exact  Jacobi  rotation  which 
transforms  H'^  to  Hm+\  in  the  commuting  diagram  of  Theorem  3.2:  jJ^H'^Jm  =  Hm+i- 

We  will  use  the  approximation  that  relgapx,  is  the  same  for  all  Hm^  even  though  it 
changes  slightly.  This  contributes  an  0{s^)  term  to  the  overall  bound  (which  we  ignore), 
but  could  be  accounted  for  using  the  bounds  of  Theorem  3.2. 

Initially,  we  will  compute  error  bounds  for  the  columns  of  Jq-  ■  ■  J\f^\,  ignoring  any 
rounding  errors  occurring  in  computing  their  product.  Then  we  wiU  incorporate  these 
rounding  errors. 

We  will  prove  by  induction  that  the  t-th  column  Vm^  of  Vm  =  Jm.-  •  ■  Jm-\  is  a  good 
approximation  to  the  true  i-th  eigenvector  I'rmi  of  Hm-  In  particular,  we  will  show  that  to 
first  order  in  e 

..  ,,    ^  (n-l)i/2(n-<o/  + M-(182(2n- 4)1/2 +  73)£)^ 

\\vTx  -  i'OtlU  <  ; 

relgapx. 

The  basis  of  the  induction  is  as  follows.  Vf^i  =  /  is  the  eigenvector  matrix  for  H^,  which 
is  considered  diagonal  since  it  satisfies  the  stopping  criterion.  Thus  the  norm  error  in  v\f^ 
follows  from  plugging  the  stopping  criterion  into  Theorem  2.7: 

II  II    ^  [n-iyl^-n-tol-R 

re/papA, 
For  the  induction  step  we  assume  that 

,,  ,,    ^  (n-l)i/2(„.  to/ +  (M-m- 1) -(182(271 -4)1/2 +  73)£)k 

\\VT,m+l,i  -  Vm+l,i\\2  < ; 

and  try  to  extend  to  m.  Consider  the  commuting  diagram  of  Theorem  3.2.  Accordingly, 
the  errors  in  V'^  =  JmKn+i  considered  ais  eigenvectors  of  H'^  are  just  the  errors  in  Vm+i 
premultiplied  by  Jm-  This  does  not  increase  them  in  2-norm,  since  Jm  is  orthogonal.  Now  we 
change  H^  to  Hm-  This  increases  the  norm  error  in  Vmi  by  an  amount  bounded  by  plugging 
the  bound  for  \\6Am\\2  into  Theorem  2.7:  (n  -  l)^/^Kil82i2n- 4y^'^ +  73)£/relgapx,-  This 
proves  the  induction  step. 
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Finally,  consider  the  errors  from  accumulating  the  product  of  slightly  wrong  values  of 
Jm  in  floating  point  arithmetic.  From  the  proof  of  Theorem  3.2,  we  see  the  relative  errors 
in  the  entries  of  Jm  are  at  most  305,  and  from  the  usual  error  analysis  of  a  product  of  2 
by  2  rotations,  we  get  32i/2.Uc  <  46-Uf  for  the  norm  error  in  the  product  of  M  rotations. 
This  completes  the  proof  of  bound  (3.12).      I 

Now  we  consider  the  errors  in  the  individual  components  of  the  computed  eigenvectors 
U'Ti(j)  -  i"i(j)l-  From  Proposition  2.12  we  see  that  we  can  hope  to  bound  this  quantity  by 
0[e]kl\{j)l  min(re/gapA, ,2"'/^).  where 

,0,.."=m,„,(ii)"^(^)'",  ,3.13, 

is  a  modified  upper  bound  for  the  eigenvector  component  i\(j)  as  in  Proposition  2.11. 
In  other  words,  we  may  have  high  relative  accuracy  even  in  the  tiny  components  of  the 
computed  eigenvectors;  this  is  the  case  in  the  example  in  the  introduction  and  at  the  end 
of  subsection  2.1.  Our  proof  of  this  fact  will  not  be  as  satisfactory  as  the  previous  result, 
because  it  will  contain  a  "pivot  growth"  factor  which  probably  grows  at  most  linearly  in  M 
but  for  which  we  can  only  prove  an  exponenticd  bound.  In  numerical  experiments  presented 
in  section  7,  there  was  no  evidence  that  this  factor  grew  with  increasing  n  or  M . 

We  will  use  v,(j)  as  defined  in  (3.13)  for  each  /?„,,  even  though  the  values  of  A,  and  Aj 
vary  slightly  from  step  to  step.  This  error  will  contribute  an  0[e^)  term  to  the  overall  bound 
(we  are  ignoring  such  terms)  but  could  be  incorporated  using  the  bounds  of  Corollary  3.4. 

Theorem  3.14  Lei  V .  Vj,  and  R  be  as  in  Theorem  3.11.  and  v,{j)  be  as  in  (3.13).  Then 
we  can  bound  the  error  in  the  individual  eigencomponents  by 

\vTi{j)-v^U)\<p{M,n)-      .  — — -  (3.15) 

nun(re/5apA,,  2-1/2) 

Here  p{M,  n)  is  a  ''pivot  growth"  factor  which  probably  grows  at  most  linearly  in  M  and  in 
n'^^^  although  all  we  can  prove  is  an  exponential  bound. 

In  numerical  experiments  presented  in  section  7,  there  was  no  evidence  that  the  factor 
p{M,n)  in  (3.15)  grew  with  increasing  n  or  M. 

Proof.  The  proof  is  similar  to  that  of  Theorem  3.11.  One  difference  is  that  we  use 
Proposition  2.12  instead  of  Theorem  2.7  to  bound  the  errors  in  the  eigenvectors.  Another 
difference,  which  introduces  the  growth  factor  p(A/),  is  that  we  need  to  use  the  scaling  of 
the  entries  of  Jm  to  see  how  small  eigenvector  components  have  small  errors;  not  being  able 
to  use  the  orthogonality  of  Jm  introduces  p(M). 

As  in  the  proof  of  Theorem  3.11,  let  V^  =  Jm  ■  ■  -Jm-Xi  where  JmH'mJm  =  Hm+i-  Set 
V\i  =  I.  The  proof  has  three  parts.  In  the  first  part  we  will  show  the  i-th  column  of  Vq  is 
a  good  approximation  to  the  eigenvectors  of  Hq  in  the  sense  of  the  theorem.  In  the  second 
part  we  will  show  that  the  (i,j)  entry  oi  Jo-  ■  ■  Jm  is  bounded  by  a  modest  multiple  of  v,(j)- 
In  the  third  part  we  will  show  the  rounding  errors  committed  in  computing  Jo---Jm  in 
floating  point  are  sm«dl  compared  to  ^.(j). 
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For  the  first  part  of  the  proof  we  will  use  induction  to  prove  that  the  z-th  column  Vm, 
of  Vm  is  a  good  approximation  to  the  true  i-th  eigenvector  VTmi  of  Hm-  This  will  show 


|iT.(J)  -  i'o,(j)l  <  Po 


'tol  + 


K  •  i\{j , 


mm(relgap\^,'2   '/'^) 


(3.161 


where  po  is  a  constant  (part  of  the  "pivot  growth"  factor)  we  need  to  estimate.  The  base 
of  the  induction  follows  from  plugging  the  stopping  criterion  into  the  bound  of  Proposition 
2.12,  yielding 


\vTMiij)  -  l^^^,{j)\  < 


2n  -  2)1/2  -n  -tot  ■  K-  i\(j] 
min(  r€lgap\,,2~^^^) 


<  PM 


{t0l  +  £)-K-  V,(J) 


where  p\f  =  n(2n  -  2)^/^.  The  induction  step  will  assume  that 


U'T.m  +  l..(i)  -   t'm  +  l,-(i)|   <  Pm  +  l 


{t0l+£] 


Mj] 


minirelgapx,,  2~^^^) 


which  we  will  try  to  extend  to  m.  Consider  the  commuting  diagram  of  Theorem  3.2. 
Accordingly,  the  errors  in  the  columns  of  V'^  =  Jm^'m+i  considered  as  eigenvectors  of  H'^ 
are  just  the  errors  in  V'^+i  premultiplied  by  Jm\  let  em,  denote  this  error  for  the  t'-th 
column  of  Vm.  Suppose  Jm  rotates  in  rows  and  columns  k  and  /;  then  em,  is  identical  to 
VT.m+i.,  -  Vm+i.,  except  for  em,ik)  and  em,{l)-  We  may  assume  without  loss  of  generality 
that  k  <  I  and  dk  >  di  {dj.  and  dj  are  the  diagonal  entries  of  Hm)-  As  in  the  proof  of 
Theorem  3.2,  there  are  two  cases,  when  x  =  di/dk  <  x  =  (\/5  -  l)/2,  and  x  >  x. 

In  the  first  case,  x  <  z,  we  know  as  in  the  proof  of  Theorem  3.2  that  s'n,  the  sine  in  the 
rotation  Jm-,  is  bounded  in  magnitude  by  i/(l  -  x'^).   Write  \s'n\  <  Cmi^i/^k)^^'^  instead. 

1  /2 

where  Cm  is  a  modest  constant.  We  can  do  this  because  dr  ss  \r  from  Proposition  2.10. 
This  lets  us  bound 


|e,n.(fc)| 
\em,{l)\ 


< 


< 


Jr, 


VT.m  +  lA^)  -  Vm+\A^) 
VT.m  +  lJl)-  Vm  +  l.x(l) 


|t^r,m+l,.(^')  -  J^m  +  l,.(fc)|  +  |sn(t;T,m  +  l,.(0  "  ^'m  +  L.CO)! 
\s'n(VT,m  +  lAf^)  -  ^m  +  \,,i'^))\  +  U'r,m  +  1,.(0  "  l'm  +  l..(OI 


min(re/pap>,,  2~i/2) 

1/2     /A.\l/2 


^ni{^r\(^y'')^cm{^s^'^M{^x^\{tr) 

c.(it)-min((^)^^(i^)-).m.n((^A.)V^(^y/^^ 
il  +  Cm)nnni{^y'\{^y^') 


Pm  +  l(t0l  +  S)K^^^ 


m.n(.e/,ap..,2-/2)     ^  (,  ^  ,^)^„(  (  A.^/^  ^  (>.y/^) 


=       (1  +  Cm) 


Pm  +  litol  +  e)K 

min(re/5apA,,  2~i/2) 


v,ik) 
r.(/) 


(3.17) 
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Now  consider  case  2,  x  >  j.  Now  A;,  and  A;  are  reasonably  close  together.  Thus,  we 
may  bound  \s'n\  simply  by  1  in  the  derivation  (3.17).  This  leads  to  the  same  bound  with 
a  possibly  different  Cm;  we  tal<e  the  final  Cm  as  the  ma.ximum  of  these  two  values.  This 
bounds  the  error  in  the  columns  of  1^  considered  as  eigenvectors  of  H'^. 

Now  we  change  //^  to  Hm-  This  increases  the  bound  for  jiTm.lj)  -  i'm,{j)\  by  an 
amount  bounded  by  plugging  the  bound  for  ||(? .4^112  from  Theorem  3.2  into  Proposition 
2.12:  (2n  -  2)^/'^(  182(2/1  -  4)'/^  +  73)  •  5  •  k-  i',(j)/ min(re/gap.\,,  2-'/^).  This  completes  the 
induction  with 


\VTn 


<     ((1  +  c,n)pm+x  +  (2n  -  2)i/2(182(2n  -  4)'/^  +  73))  • 


=        Pr 


min(re/5apA,, 2-1/2) 


{t0l  +  £)-K-V,{j) 

m\n(relgap\^,  2~'^l'^) 
(3.18) 


Here 


=  (l  +  c^)p„  +  i+(2n-2)'/'(182(2n-4)'/'  +  73)    .    pM  =  n[2n-2)"^ 


(3.19) 


satisfies  an  exponential  error  bound,  but  it  is  clear  from  the  derivation  that  linear  growth 
is  far  more  likely  than  exponential  growth.  This  completes  the  first  part  of  the  proof. 

In  the  second  part  of-the  proof  we  will  show  that  the  (j,j)  entry  of  V^  =  Jo  ■  ■  ■  Jm  is 
bounded  by  a  modest  multiple  of  t\(j).  To  do  this  we  will  prove  by  induction  that 


\i'm.U]\  <  TmViU) 


(3.20) 


where  V^  =  [i'mi  ,••■•.  Vmn]  and  r^  is  a  constant  (part  of  the  ''pivot  growth"  factor)  we  need 
to  estimate.  The  base  of  the  induction  is  for  m  =  -1,  i.e.  the  null  product,  which  we  set 
equal  to  the  identity  matrix.  This  clearly  satisfies  (3.20)  with  r_i  =  1.  Now  we  assume 
(3.20)  is  true  for  m  —  1  and  try  to  extend  it  to  m.  Suppose  Jm  rotates  in  rows  and  columns 
k  and  /.  Postmultiplying  Vm-i  by  J^  only  changes  it  in  columns  k  and  /.  Assume  as  before 
that  k  <  I  and  x  =  di/dk  <  1.  There  are  two  cases  as  before,  i  <  i  and  x  >  x. 

First  consider  the  case  x  <  x.    We  may  bound  the  (j,k)  and  (j,l)  entries  of  V„,  as 
follows: 


IK,,(fc),t'm,,(/)]|-    =       \[Vm-lAk),Vm-uU)] 


cs 
—s'n 


< 


^        '  m 


sn 
cs 

1/2 


),min((^ 


1/2 


?         )]• 


< 


=       T, 


=        Tr 


1  -  (It)'" 

/\    \  1/2      /\,  \^/2  /\    N  1/2 

,-lil  +  Cm)[vk{j),Vl{j)] 

,[vk{j),vi(j)] 


1/2 


(3.21) 
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In  the  second  case,  x  >  x,  we  get  a  similar  bound.  Here  Xk  ~  -^z-  ^nd  we  can  simply  bound 
\s~n\  <  1.  This  yields  a  slightly  different  c^;  for  the  final  Cm  we  again  take  the  maximum  of 
the  two.  This  ends  the  second  part  of  the  proof  with 

Tm   =   (I  +  Cm)Tm-i     ,     r_i    =    l  (3.22) 

Even  though  this  only  yields  an  exponential  upper  bound  for  rv/,  it  is  clear  from  the 
derivation  that  linear  growth  is  far  more  likely  than  exponential  growth. 

In  the  third  and  final  part  of  the  proof  we  show  the  rounding  errors  in  the  (i.j)  entry 
of  the  computed  approximation  to  Vm-\  is  bounded  by  0(£)v,{j).  Let  Jm  be  the  ac;ual 
rotation  which  only  approximates  Jm-  From  the  proof  of  Theorem  3.2,  we  have  cs  = 
cs{l  +  £cs)  with  Iscsl  <  225  and  sn  =  5n(  !+£,„)  with  |£,„|  <  30£-.  Let  Vm  =  IH^'m-x  *  Jm) 
be  the  actually  computed  eigenvector  matrix  after  the  m-th  Jacobi  rotation.  The  final 
computed  eigenvector  matrix  is  V  =  V',\/_i.  We  will  use  induction  to  prove 

I'vm.Aj)  -  VmAJ)\  <   \mSV,U)  (3.23) 

where  V'^  =  [I'mi.  •  •  • ,  u^n]  and  \m  is  a  constant  (part  of  the  pivot  growth"  factor)  we 
need  to  estimate.  The  basis  is  again  for  m  =  —  1  when  V^i  =  VLj  =  /  and  \_i  =  0.  Now 
we  assume  (3.23)  is  true  for  ttj  —  1  and  try  to  extend  it  to  m.  As  before,  we  assume  J^ 
rotates  in  rows  and  columns  k  and  /  with  k  <  /  and  x  =  dk/di  <  1.  Write  €„,  =  Vmi  -  Vmi- 
The  {j,k)  and  (j, /)  entries  of  \'m  are 

[Vmj{k)        ,        Um,;(0] 

=       [C„^_i,j(fc)c3(l  +  £cs)(l  +  fl)(l  +  ^2)  -  'Vm-i,:{^)sn{l  +  £,n){l  +  f3)(l  +  ^4), 

'Vm--ij(k)s-n(l+  f,n)(l  +  £5)(1  +  Ss)  +  'Vm-l.j{l)cs(i  +  £c,)(l  +  f7)(l  +  £»)] 
=       [Vm-l.j(k)CS  -  Vm-l.j{l)sn,  Vm-\,j{k)sn  +  V'„,_i,j(/)C5]  + 

[2A£sCSVm-\,j{k)  +  32£iQSnVm-l,jil),32£iiSnVm-\j{k)  +  24Si2CSVm-l.j(l)]  + 

[(1  +  24£9)c3em-i,:(k)  +  (1  +  32£ioKne„_i,;(/), 
(1  +  32£ii)5ne„_i,j(fc)  +  (1  +  24ei2)c5e„-i,;(/)] 
=     [i;„,^(fc),r^,j(/)]  +  /i+/2 

so  [em,](k),em.]ih]  =  h  +  h- 

As  before,  there  are  two  ca^es,  i  <  i  and  x  >  x.  Consider  the  case  x  <  x.  Using 
\sn\  <  CmiXi/Xk^^'^,  \cs\  <  1,  and  |i;m-i,.(j)|  <  7-„,_iv,(j),  we  get 

1/2  .,     X  1/2       /  V     \  1/2 


X 


/A  \^/^    /a   \  /A 

l/il     <     £r„-iK^/'[(24  +  32c^)min((^^j       ,(^j       ),  (24  +  32c„)  min(  (^^ 

=      £7^-1(24  +  32Cm)[vkU),Vk{l)] 

and 

Ihl  <  Xm-l(l  +  Cm)e[VkU)rVk{l)] 

Taken  together,  we  get 

\m  =  (l  +  c,„)Xm-i +£T„-i(24  +  32c„)   ,    x-i=0  (3.24) 
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In  the  second  case,  x  >  x.  we  get  a  similar  bound  with  a  possibly  different  c^.  Again,  we 
take  the  maximum  of  the  two.  This  completes  the  third  part  of  the  proof. 
Finally,  combining  (.3.23)  and  (3.16)  we  get 

lM.)-rr,.(;)|<(Po+X.u-.^    (fo/  +  ,)-^-r.(;) 


min(  relgap\,,  2~'^l'^) 
proving  the  theorem  with  p(M,n)  =  po  +  \.\/-i.       I 
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4      One-sided  Jacobi 

In  this  section  we  prove  that  one-sided  Jacobi  in  floating  point  arithmetic  applied  to  a 
general  matrix  computes  the  singular  values  and  singular  vectors  with  the  error  bounds  of 
section  2.  Here  we  present  our  algorithm;  the  model  of  arithmetic  was  presented  in  section 
3.  In  subsection  4.1  we  derive  error  bounds  for  the  computed  singular  values.  In  subsection 
4.2  we  derive  error  bounds  for  the  computed  singular  vectors.  In  subsection  4.3,  we  present 
two  algorithms  for  the  symmetric  positive  definite  eijenproblem  H .  both  of  which  involve 
applying  one-sided  Jacobi  to  the  Cholesky  factor  L  (or  L)  of  H.  The  second  of  these 
algorithms  annot  compute  eigenvectors  quite  as  accurately  as  the  first,  but  may  be  much 
faster  than  either  the  first  adgorithm  or  two-sided  Jacobi. 

Let  Go  =  BqDq  be  the  initial  matrix,  and  Gm  =  BmDm,  where  Gm  is  obtained  from 
G'm-i  by  applying  a  single  Jacobi  rotation.  Here  Dm  is  diagonal  and  Bm  has  columns  of 
unit  norm.  All  the  error  bounds  in  this  section  contain  the  factor  maxm'^(5m),  whereas 
the  perturbation  bounds  in  section  2  are  proportional  to  k{Bq).  Therefore,  as  in  section  3, 
our  claim  that  Jacobi  computes  the  SV'D  as  accurately  as  predicted  in  section  2  depends 
on  the  ratio  ma^x-m '^{Bm)l i^{Bq)  being  modest.  In  exact  arithmetic,  one-sided  Jacobi  on 
G  =  BD  is  identical  to  two-sided  Jacobi  on  H  =  G^G  =  DB^ BD  =  DAD,  so  the  question 
of  the  growth  of  ^(5^  )  =  KiAm)^^^  is  essentially  identical  to  the  question  of  the  growth  of 
K{Am)  in  the  case  of  two-sided  Jacobi. 

The  essential  difference  between  our  algorithm  and  standard  one-sided  Jacobi  is  the 
stopping  criterion:  according  to  Theorem  2.17,  we  must  stop  when  all  n^jl{H„Hjj)^l'^ 
are  small  {H  =  G^G),  not  just  when  Htj/  maxkilNkil  is  small.  This  stopping  criterion 
has  been  suggested  before  [16,5,3],  but  without  our  explanation  of  its  benefits.  Otherwise, 
our  algorithm  is  based  on  the  standard  one  introduced  by  Rutishauser  [12].  We  have 
chosen  a  simple  version  of  the  algorithm,  omitting  enhancements  like  delayed  updates  of 
the  diagonals  and  fast  rotations,  to  make  the  error  analysis  clearer  (an  error  analysis  of 
these  enhancements  is  future  work). 

Algorithm  4.1  One-sided  Jacobi  for  the  general  singular  value  problem,  tol  is  a  use  defined 
stopping  criterion.  The  matrix  V  whose  columns  are  the  computed  right  singular  vectors 
initially  contains  the  identity. 


repeat 

for  all  pairs  i  <  j 

/*  compute 

«  =  ELi  Gl 
b  =  TIU  Gl, 


a    c 

c    b 


=  the  {i,j)  submatrix  of  G^G  */ 


c  =  TX=\  Gk,  *  Gkj 

/*  compute  the  Jacobi  rotation  which  diagonalizes 

C  =  {b-a)/{2c) 


t  =  sign{o/(\c\  +  yr+c^) 
cs  =  i/>/rT^ 


a    c 

c    b 


V 
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sn  =  cs  *  t 

/*  update  columns  i  and  j  of  G  */ 
for  k  =  I  to  n 
tmp  =  Gkx 

Gki  =  cs  ♦  tmp  -  sn  *  Gkj 
Gkj  =  sn  ♦  tmp  +  cs  *  Gkj 
endfor 

/*  update  the  matrix  \'  of  right  singular  vectors  */ 
for  k  —  I  to  n 
tmp=  Vkx 

Vkx  =  C5  *  tmp  -  sn*  Vkj 
Vkj  =  sn  *  tmp  -\-  cs  *  V'kj 
endfor 
endfor 
until  convergence  (all  \c\l\/ab  <  tol) 

/*  the  computed  singular  values  are  the  norms  of  the  columns  of  the  final  G  */ 
/*  the  computed  left  singular  vectors  are  the  normalized  columns  of  the  final  G  */ 

4.1      Error  Bounds  for  Singular  Values  Computed  by  One-sided  Jacobi 

The  next  theorem  and  its  corollary  justify  our  accuracy  claims  for  singular  values  computed 
by  one-sided  Jacobi. 

Theorem  4.1  Let  Gm  ^  '/^e  sequence  of  matrices  generated  by  the  one-sided  Jacobi  algo- 
rithm in  finite  precision  arithmetic  with  precision  e;  that  is  Gm+i  ^s  obtained  from  Gm  by 
applying  a  single  Jacobi  rotation.    Then  the  following  diagram  commutes. 

G      floating     /~T 

+SG. 


The  top  arrow  indicates  that  Gm+i  ^s  obtained  from  Gm  by  applying  one  Jacobi  rotation 
in  floating  point  arithmetic.  The  diagonal  arrow  indicates  that  Gm+i  is  obtained  from 
G'm  by  applying  one  Givens  rotation  in  exact  arithmetic:  thus  Gm+i  omd  G'^  have  identical 
singular  values  and  left  singular  vectors.  The  vertical  arrow  indicates  that  G'm  =  Gm  +  ^Gm- 
6Gm  is  bounded  as  follows.  Write  6Gm  —  ^BmDm,  where  Dm  's  diagonal  such  that  Bm  >" 
Gm  =  BmDm  hos  Unit  columns.   Then 

WSBmh  <  72e  (4.2) 

In  other  words,   one  step  of  Jacobi  satisfies  the  assumptions  needed  the  error  bounds  of 
section  2. 

Corollary  4.3  Assume  Algorithm  4-1  converges,  and  that  Gm  is  the  final  matrix  which 
satisfies  the  stopping  criterion.  For  0  <  m  <  M  write  Gm  =  BmDm  with  Dm  diagonal  and 
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Bm   with  unit  columns.    Let  a^   he  the  j-th  singular  value  of  Gq  and  (t'   the  j-th  computed 
singular  value.    Then  to  first  order  in  s  the  following  error  bound  holds: 

I  _    _  n-'  I 

— <  {72s  ■  M  +  n'^s  +  n-tol)  ■    max   K(Bk)  +  n£  (4.4) 

<7j  ~  0<k<.\{ 

Remark.  In  numerical  experiments  presented  in  section  7,  there  was  no  evidence  that  the 
error  bound  in  (4.4)  grew  with  increasing  n  or  M. 

Proof  of  Corollary  4.3.   Bound  (4.4)  follows  by  substituting  the  bound  (4.2)  and 
the  stopping  criterion  into  Theorem  2.17.   The  n'' r  term  comes  from  the  fact  that  c/y/ab 
in  the  stopping  criterion  may  be  underestimated  by  as  much  as  ne.  The  trailing  ns  comes 
from  computing  the  norms  of  the  columns  of  the  final  G  matrix.       I 
Remark.  A  similar  bound  can  be  obtained  based  on  the  error  bound  in  Proposition  2.19. 

Proof  of  Theorem  4.1.  The  proof  of  the  commuting  diagram  is  a  tedious  computa- 
tion. Let  aj,  bj  and  cj  be  the  true  values  ofYlk^li^  Hi  ^L  ^"d  "^Zk^kiG^j.  Then  (in 
the  notation  of  the  proof  of  Theorem  3.2)  write  a-p  =  d^,  bj-  =  d^  and  cj  —  :d,dj.  We 
may  assume  without  loss  of  generality  that  aj  >  bj  and  cj  >  0.   By  positive  definiteness 

0  <  z  <  :  =  {K^iBm)  -  l)/{K^{Bm)  +  1)  <  1.   Let  X  =  dj/d,  <  1.  We  consider  two  cases, 

1  <  i  =  (v/5  -  l)/2  %  .62,  and  X  >  X. 

First  consider  x  <  x.  Systematic  application  of  formulas  (3.1)  shows  that 

a  =  a j(  1  +  fa  )     where     |fa|  <  "f 

b  =  bjil  +  Bb)     where     \sb\  <  ne 

c  =  CT  +  ec\/ajbj     where     \Ec\  <  ne 

Let  C3  =  (l  +  <2)-i/2  and  s'n  =  <(1  +  t2)-i/2.  Then  from  (3.1 )  again  we  get  sn  =  {l  +  £,r^)sn 
and  C5  =  ( 1  +  Sc3)cs  where  |5,„|  <  4£  and  \£c,\  <  Se.  c's  and  s'n  define  the  Givens  rotation 
cs  s'n 
—s'n  c's 
(1  +  0{£))x{z  +  n£)/(l  -  x^)  and  so  |5n|  <  {I  +  0(£))x{z  +  n£)/(l  -  x^). 

Let  G'^,  and  6"^^  be  the  new  values  for  these  entries  computed  by  the  algorithm.  Using 
the  bounds  |c5|  <  1  and  \s'n\  <  (1  -I-  0(£))x{z  +  n£)/{l  -  i^),  we  estimate 

G'fc,     =  //(cs  *  Gfc,  -  sn  ♦  Gfcj) 

=  (1  +  £i)(l  +  £2)csGk,  -  (1  +  £3)(1  +  S4)snGkj 

=  (1  +  £i)(l  +  £2)(1  +  £cs)csGk,  -  (1  +  £3)(1  +  f4)(l  +  e,n)snGk: 

=  c'sGkt  -  s'nGkj  +  Eki  (4.5) 


'^TTl      — 


which  take  G^  to  Gm+x'-   G^Jm  =  Gm+i-    Also,  we  can  show  t  < 


and 


G'kj     =  fl{sn*Gk,  +  cs*Gkj) 

=  (1  +  £5)(1  +  £6)snGk,  +  (1  +  S7){1  +  Ss)csGki 

=  (1  +  £5)(1  +  £6)(1  +  e,r,)s'nGk,  +  (1  +  £7)(1  +  e8)(l  +  £c.)csGk: 

=  s'nGkt  +  csGk:  +  Ekj  (4.6) 
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where 


1^.112     <     5£||G.,||2  +  6£||G.,||2  <  Usd, 


6{z+ns) 
[here  G.,  refers  to  the  i-th  column  of  G.  etc.).  Thus 


+  5)dj 


cs       sn 

-  s'n     cs 


G\    G',  ]     =     [  G.,    G., 

G'.,     G.,  J  +  [  £..     Ej 
G..    G.,  ]  +  [  F.     F.,  ]) 


+     F.    £.j 


C3  —sn 

s'n  cs 

cs  s'n 

-  s'n  cs 


cs       sn 
—  s'n     cs 


where 


\F.,\\2     <     ||£..||2  +  l|F.,|b<£(^Y^^  +  16K 
\F.A\2     <     \s'n\-\\E.,h  +  \\Ej2<s(^^^^^  +  5)d, 


Thus 


l-Smlb   < 


^  +  ^<,,?Ei±2£)+21,  ,4.7) 

a,  dj  1  -  x^ 

Now  consider  the  case  x  >  x.  The  analysis  differs  from  the  previous  one  only  in  the  fact 
that  sn  is  no  longer  small.  Using  the  bounds  \s'n\  <  1,  |c5|  <  1  in  (4.5)  and  (4.6)  yields 

||£..||2     <    5£||G.,||2  +  6e||G.j||2  <  lied. 

||£.,1|2      <      6£||G.,||2  +  5£||G.;||2<ll£d. 

whence 

IIF.II2     <     \\E.,\\2  +  \\E.j\\2<22ed, 

\\F.,\\2       <       |i£.,||2  +  ||£.,||2<22£(i,/x 


and 


|5.IU<^  +  ^<44^/x 


d, 


d: 


(4.8) 


Since  x  satisfies  1/(1  —  i^)  =  l/i  <  1.62,  we  see  from  (4.7)  and  (4.8)  that  in  both  cases 


Bmh  <  72£ 


proving  the  theorem.      I 
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4.2      Error  Bounds  for  Singular  Vectors  Computed  by  One-sided  Jacobi 

The  next  two  theorems  justify  our  accuracy  claims  for  singular  vectors  computed  bv  one- 
sided Jacobi. 

Theorem  4.9  Let  V  =  [I'l,  ■  •  •,  r,,]  be  the  matrix  of  unit  right  singular  vectors  and  L'  = 
[ui.-  ■  -.Un]  be  the  matrix  of  unit  left  singular  vectors  computed  by  Algorithm  4.I  in  finite 
precision  arithmetic  with  precision  e.  Let  X'r  =  [i-'ji,  •  •  • ,  UTn]  and  Uj  =  [uji- ■  •  •  ■  "Tn]  be 
the  matrices  of  true  unit  right  and  left  singular  vectors,  respectively.  Let  n  =  max^  K(Bm) 
be  the  largest  K(Bm)  of  any  iterate.  Then  the  error  in  the  computed  singular  vectors  ■■ 
bounded  in  norm  by 

,,,  ,,     .,  „,       (n-  .0)^^^ -R-ilQM -s +  n-tol +  n^  ■£) 

max(||ur,  -  u,||2,|ri'T.  -  t'.l  2)  <  ; +  (9-U  +  n  +  l)e 

relgapa^ 

(4.10) 

Remark.  In  numerical  experiments  presented  in  section  7,  there  was  no  evidence  that  the 
error  bound  in  (4.10)  grew  with  increasing  n  or  .\/. 

Proof.  The  proof  is  similar  to  that  of  Theorem  3.11.  Let  G'o,  • .  •  ,<j.vf  be  the  sequence 
of  matrices  generated  by  the  Jacobi  algorithm,  where  Gm  satisfies  the  stopping  criterion. 
Let  Jrn  be  the  exact  Givens  rotation  which  transforms  G^  to  Gm+\  in  the  commuting 
diagram  of  Theorem  4.1:  G'^Jm.  =  Gm+i- 

We  will  use  the  approximation  that  relgap„^  is  the  same  for  all  Gm,  e\cu  though  it 
changes  slightly.  This  contributes  an  0(£^)  term  (which  we  ignore),  but  could  be  accounted 
for  using  the  bounds  of  Theorem  4.1. 

First  we  consider  the  left  singular  vectors.  In  exact  arithmetic,  these  remain  unchanged 
throughout  the  computation  since  all  rotations  are  applied  on  the  right.  Thus,  we  need 
only  plug  the  bounds  for  the  stopping  criterion  and  each  ||(55rn||2  from  Theorem  4.1  into 
Theorem  2.21  to  get 

wr.  -  w.   2  <  ^ -, +  (n  +  1  £ 

Telgap^i 

as  claimed  (the  (n  +  1)£  term  comes  from  normalizing  the  i-th  column  of  G.v/  at  the  end  of 
the  computation).* 

Now  we  consider  the  right  singular  vectors.  First  we  will  compute  error  bounds  for  the 
columns  oi  Jq-  ■  •  J\i-\  ignoring  any  rounding  errors  occurring  in  computing  their  product. 
Then  we  will  incorporate  these  rounding  errors. 

We  wiU  prove  by  induction  that  the  i-th  column  Vmt  of  Vi„  =  7^  •  •  Jm-i  is  a  good 
approximation  to  the  true  i-th  right  singular  vector  ut-^,  of  Gm-  In  particular,  we  wll  show 
that  to  first  order  in  e 

{n- .5V^^-K-{72M -e  +  n-tol  +  n^  ■£) 
\\VTi  -  I'D.    2  <  ; 

The  basis  of  the  induction  is  as  follows.  V\f  =  /  the  the  singular  vector  matrix  for  Gm, 
which  is  considered  to  have  orthogonal  columns  since  it  passes  the  stopping  criterion.  Thus 
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the  norm  error  in  v\f,  follows  from  plugging  the  stopping  criterion  n  •  tol  (increased  by  n^s 
since  n  •  tol  may  be  underestimated  by  this  amount)  into  Theorem  2.21: 

„  „         {n-  .5)^^'^  -R-in-tol  +  n^  -s) 

For  the  induction  step  we  assume  that 

(n  -  .5)i/2-«-(72(.\/  -  m  -  I)  ■  E  +  n  ■  tol  +  n^  ■  s) 


h'T.m  +  l.i  —   i-'m  +  l.ilb   £ 


relgapa, 


and  try  to  extend  to  m.  Consider  the  commuting  diagram  of  Theorem  4.1.  Accordingly, 
the  errors  in  I'm  =  Jm^m  +  \  considered  as  right  singular  vectors  of  G'^  are  just  the  errors 
in  V-m+i  premultiplied  by  Jm-  This  does  not  change  their  norm,  since  Jm  is  orthogonal. 
Now  we  change  G^  to  Gm-  This  increases  the  norm  error  in  Vm,  by  an  amount  bounded  by 
plugging  the  bound  for  ||<^5m||2  into  Theorem  2.21:  12e{n  -  .Sy^^R/relgap^,.  This  proves 
the  induction  step. 

Finally,  consider  the  errors  from  accumulating  the  product  of  slightly  wrong  values  of 
Jm  in  floating  point  arithmetic.  From  the  proof  of  Theorem  4.1,  we  see  the  relative  errors 
in  the  entries  of  Jm  are  at  most  -Is,  and  from  the  usual  error  analysis  of  a  product  of  2  by 
2  rotations,  we  get  6\/2M^  <  9Me  for  the  norm  error  in  the  product  of  M  rotations.  This 
completes  the  proof  of  the  theorem.       I 

Now  we  consider  the  errors  in  the  individual  components  of  the  computed  right  singular 
vectors  \vti{j)  -  ^\[j)\-  From  Proposition  2.26  we  see  that  we  can  hope  to  bound  this 
quantity  by  0{s)Rv,(j)/relgapa,,  where 

v,(j)  =  K^mm(—,^)  (4.11) 

aj    a, 

is  a  modified  upper  bound  for  the  right  singular  vector  component  v,{j)  as  in  Proposition 
2.25.  In  other  words,  we  may  have  high  relative  accuracy  even  in  the  tiny  components 
of  the  computed  right  singular  vectors.  Our  proof  of  this  fact  will  not  be  as  satisfactory 
as  the  previous  result,  because  it  will  contciin  a  "pivot  growth"  factor  which  probably 
grows  at  most  linearly  in  M  but  for  which  we  can  only  prove  an  exponential  bound.  In 
numerical  experiments  presented  in  section  7,  there  was  no  evidence  that  this  factor  grew 
with  increasing  n  or  M . 

We  will  use  v,(j)  as  defined  in  (4.11)  for  each  Gm-,  even  though  the  values  of  ct,  and 
a  J  vary  slightly  from  step  to  step.  This  error  will  contribute  an  O(s^)  term  to  the  overall 
bound  (which  we  ignore)  but  could  be  incorporated  using  the  bounds  of  Corollary  4.3. 

Theorem  4.12  Let  V,  Vj  and  k  be  as  in  Theorem  4-9,  and  v,(j)  be  as  in  (4.11).  Then 
we  can  bound  the  error  in  the  individual  components  of  v,  by 

MJ) -  v.(;)l  <  .(M, .) .  i^oUe).R^.v,u)  ^^^3^ 

relgap^, 

Proof.  The  proof  is  nearly  identical  to  that  of  Theorem  3.14;  we  just  outline  the 
differences  here.    Let  Jm  be  the  exact  Givens  rotation  whic  transforms  G'm  to  Gm+i  in 
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the  commuting  diagram  of  Tlieorem  4.1:  G'^Jm  =  G^r.+i-  Let  V'^  =  J^---Jvr-i  and 
V'^  =  Jo---Jm-  In  the  first  part  of  the  proof  we  show  the  columns  of  V'o  have  small 
componentwise  errors  in  the  sense  of  the  theorem.  In  the  second  part  we  will  show  each 
(i.j)  entry  of  each  V'^  is  bounded  by  a  modest  multiple  of  v,(j).  In  the  third  part  we  show 
the  rounding  errors  committed  in  computing  I'm  in  floating  point  are  componentwise  small 
compared  to  f,(j). 

For  the  first  part,  the  same  induction  argument  as  in  Theorem  3.14  leads  to  the  bound 

\VTmAj)~   VmAj)\   <  Pm- 


relgapa, 


where 


U4^/n^^^ 


P'm    =   (  1  +  Cm  )P'm  +  l  +  ' '     p'sf   =  ^uVrT^ 

Here  c'^  is  a  small  constant  as  in  Theorem  3.14.  We  use  Proposition  2.26  and  Theorem  4.1 
in  place  of  Proposition  2.12  and  Theorem  3.2  which  were  used  in  Theorem  3.14. 

For  the  second  part,  let  IVn   =   [i'mi i'mn]-    The  same  induction  argument  as  in 

Theorem  3.14  yields 

|f'm.(j)l<<l=.(j) 

where 

^m    =  (1  +  C„,)r„_i     ,      T_i    =    1 

For  the  third  part,  let  V'^  =  /'(^'m-i  *  Jm)  be  the  actually  computed  singular  vector 

matrix  after  the  m-th  rotation.  Write  V'm  =  [vmi>  ■  ■  -,  Vmn]-  The  same  induction  argument 
as  in  Theorem  3.14  yields 

\Vm,U)-  Vm,U)\  <  x'mfMJ) 

with 

Altogether,  we  get 

I       ,  ■\         I -w  ^  (   '  ^    '        ,(<o^+  ne)-K^  -v.ij) 

\vT,{j)-  v,{j)\  <  (Po  +  Xm-i) ; 

relgap^, 

proving  the  theorem  with  q{M,  n)  =  Pq  +  Xa/-i-       ' 

4.3     Using  Cholesky  Followed  by  One-sided  Jacobi  for  the  Symmetric  Pos- 
itive Definite  Eigenproblem 

In  this  subsection  we  consider  two  algorithms  for  the  symmetric  positive  definite  eigenprob- 
lem H ,  both  based  on  performing  Cholesky  on  H,  and  using  one-sided  Jacobi  to  compute 
the  SVD  of  L.  The  first  algorithm  (Algorithm  4.2)  does  one-sided  Jacobi  on  L^ ,  returning 
its  right  singular  vectors  as  the  eigenvectors  of  H  and  the  squares  of  its  singular  values  as  the 
eigenvalues  of  H .  The  second  algorithm  (.Algorithm  4.4),  originally  proposed  in  [16],  does 
Cholesky  with  complete  pivoting,  and  then  one-sided  Jacobi  on  L,  returning  its  left  singular 
vectors  as  the  eigenvectors  of  H  and  the  squares  of  its  singulcir  values  as  the  eigenvalues 
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of  H .  The  second  algorithm,  which  we  call  accelerated  one-sided  Jacobi,  is  less  accurate 
than  the  first  because  it  will  not  always  compute  tiny  eigenvector  components  with  the 
accuracy  of  Theorem  3.14.  although  it  does  compute  the  eigenvalues  as  accurately,  and  the 
eigenvectors  with  the  same  norm  error  bound.  However,  it  can  be  several  times  faster  than 
either  the  first  algorithm  or  two-sided  Jacobi.  In  fact,  the  larger  the  range  of  numbers  on 
the  diagonal  of  D.  the  faster  the  second  algorithm  will  converge.  This  means  that  the  more 
the  guaranteed  accuracy  of  the  algorithm  exceeds  that  of  QR  (or  any  tridiagonalization 
based  algorithm),  the  faster  it  converges. 

Algorithm  4.2  One-sided  Jacobi  method  for  the  symmetric  positive  definite  eigenproblem 
H. 

1.  Form  the  Cholesky  factor  L  of  H :  H  —  LL^ . 

2.  Compute  the  singular  values  a,  and  right  singular  vectors  f,   of  L      using  one-sided 
Jacobi. 

3.  The  eigenvalues  A,  of  H  are  A,  =  a].   The  eigenvectors  of  H  are  v,. 

We  show  this  method  is  as  accurate  as  using  two-sided  Jacobi  directly  on  H .  The  proof 
involves  a  new  error  analysis  of  Cholesky  decomposition,  so  we  begin  by  restating  Cholesky 's 
algorithm  in  order  to  establish  notation  for  our  error  analysis: 

Algorithm  4.3  Cholesky  decomposition  H  =  LL^  for  an  n  by  n  symmetric  positive  definite 
matrix  H . 

for  I  =  1  to  n 

for  j  =  i  -\-  I  to  n 

Tj,  =  (Hj,  -  XIfc=i  Tjk^xk}/L„ 
endfor 
endfor 

Lemma  4.14  Let  L  be  the  Cholesky  factor  of  H  computed  using  Algorithm  4-3  in  finite 
precision  arithmetic  with  precisions.   Then  LL^  =  H  +  E  where  |£,j|  <  {n  +  5)e{HiiHjj)^^'^. 

Proof.    Applying  rules  (3.1)  for  floating  point  arithmetic  yields 

Lu  =  (1  +  £i)((l  +  £2)^..  -  E  ^^k(^  +  i^k+2))'^^ 

k=i 

Since  ^).'l\  Lji^  <  H„  to  within  a  small  relative  error,  we  may  write 

k=i 
as  desired.  Next  we  have 

1-1 

L,j  =  (l  +  5i)((l+£2)^j.-  ^£jfcI,Ul  +  2fi+2) 

k=i 
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■  1-1 

-<:=1 

inequality,  we  can  write 


since  \  ^[Ji  LjkL,k\  <  (//'„//'jj  )^''^  to  within  a  small  relative  error  by  the  Cauchy-Schwartz 


proving  the  result.       I 

Theorem  4.15  Let  L  be  the  Cholesky  factor  of  H  =  DAD  computed  in  floating  point 
arithmetic  using  Algorithm  4-3.  Let  a,  and  r^,,  be  the  exact  singular  values  and  right  singular 
vectors  of  L  .  and  A,  and  Vfj,  be  the  eigenvalues  and  eigenvectors  of  H .  Let  i\(j)  be  as  in 
Proposition  2.11.    Then 

lA.  -  a^\ 


A, 


<  (n^  +  5n)  •£ • k{A) 


relgapx, 

I      ,■^  ,  .,,^  {n^  +  5n)i2n-2)'^^-s-K{A)-v,(j)   ,   ^,  ,^ 

\vUj)-  VH:ij)\  <  -— -; —-77- -^  +  0(£') 

min(relgapx^.2    ''^) 

Proof.  Plug  the  bound  of  Lemma  4.14  into  Theorem  2.3.  Theorem  2.7  and  Proposition 
2.12.       I 

Theorem  4.15  implies  that  the  errors  introduced  by  Cholesky  are  as  small  as  those 
introduced  by  two-sided  Jacobi.  Write  H  =  DAD  and  La  =  D~^L.  Since  ||A-  LaL'^\\2  < 
(n^  +  5n)5,  k(A)  »  {k(La))^  (unless  both  are  very  large).  Since  the  columns  of  ij  have 
nearly  unit  norm,  the  accuracy  of  one-sided  Jacobi  applied  to  L^  is  governed  by  k(Z^). 
Thus,  Cholesky  followed  by  one-sided  Jacobi  results  in  a  problem  whose  condition  number 
k(La)  is  approximately  the  square  root  of  the  condition  number  of  the  origincd  problem 
k(A).  Corollary  4.3  and  Theorems  4.9  and  4.12  guarantee  that  the  computed  eigenvalues 
and  eigenvectors  are  accurate.  In  exact  arithmetic  one-sided  Jacobi  on  L^  is  the  same  as 
two-sided  Jacobi  on  DAD  =  H  =  LL^  =  D{LaL^)D,  so  the  question  of  how  much  k(La) 
can  grow  during  subsequent  Jacobi  rotations  is  essentially  identical  to  the  question  of  the 
growth  of  K(Am)  during  two-sided  Jacobi. 

Here  is  the  second  algorithm: 

Algorithm  4.4  Accelerated  one-sided  Jacobi  method  for  the  symmetric  positive  definite 
eigenproblem  H . 

1.  Form  the  Cholesky  factor  L  of  H  using  complete  pivoting.   Then  there  is  a  permutation 
matrix  P  such  that  P^ H P  =  LL'^. 

2.  Compute  the  singular  values  a,  and  left  singular  vectors  u,  of  L  using  one-sided  Jacobi. 

3.  The  eigenvalues  A,  of  H  are  A,  =  cr^.   The  eigenvectors  of  H  are  Pu,. 

Even  if  we  did  not  do  complete  pivoting,  Theorem  4.15  would  guarantee  that  the  squares 
of  the  true  singular  values  of  L  would  be  accurate  eigenvalues  of  H ,  and  that  the  true  left 
singular  vectors  of  L  would  be  accurate  eigenvectors  of  P^ HP.    Since  we  are  computing 
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left  singular  vectors  of  I,  Theorem  4.12  does  not  apply,  but  from  Corollary  4.3  we  know 
the  computed  eigenvalues  are  accurate,  and  from  Theorem  4.9  we  know  the  computed 
eigenvectors  are  accurate  in  a  norm  sense.  Numerical  experiments  in  section  7  below  bear 
out  the  fact  that  tiny  eigenvector  components  may  not  always  be  computed  as  accurately 
by  Algorithm  4.4  as  Algorithm  4.2. 

The  advantage  of  complete  pivoting  is  accelerated  convergence.  This  is  because  the 
algorithm  is  in  principle  doing  two-sided  Jacobi  on  L^ L,  so  writing  L^ L  =  H'  =  D'A'D'. 
it  is  k(A')  which  the  algorithm  must  drive  to  1.  not  n{A).  The  wider  the  range  of  numbers 
on  the  diagonal  of  D.  the  smaller  k{A')  will  be.  We  discuss  this  in  more  detail  in  section  6, 
and  content  ourselves  here  with  a  small  example.  Let 


A  = 


1 
a 


D  =  diag  (l,c/)     and       H  = 


1      da 
da     d^ 


where  0  <  a  <  1,  and  0  <  d  <  1  (any  2  by  2  symmetric  positive  definite  H  can  be  scaled 
and  permuted  to  be  in  this  form).  Here  k(A)  =  (l  +  a)/(l-a),  which  can  be  made  as  large 
as  desired  by  choosing  a  near  1.  The  matrix  is  already  ordered  for  complete  pivoting,  and 
the  Cholesky  factor  is 

1  0 

da    (f(l-a2^i/2 


L  = 


so  L'^L  =  H'  =  D'A'D'  where 


A'  = 


1  (l  +  dV)-l/2 


and  k(A')  <  (3  +  ^/(b))l2  ss  2.62  independent  of  H .  In  section  6  we  will  prove  in  general 
that  k{A')  is  bounded  independent  of  H . 
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5      Bisection  and  Inverse  Iteration 

Here  we  show  bisection  and  inverse  iteration  applied  to  the  symmetric  positive  definite 
matrix  H  =  DAD  can  compute  the  eigenvalues  and  eigenvectors  within  the  accuracy  bounds 
section  of  2.  Let  inertia(//  )  denote  the  triple  (n.  z,p)  of  the  number  n  of  negative  eigenvalues 
of  H .  the  number  ;  of  zero  eigenvalues  of  If.  and  the  number  p  of  positive  eigenvalues  of 
H .  These  results  are  extensions  of  Algorithms  3  asid  5  in  [2]. 

Algorithm  5.1   Stably  compisiing  the  inertia  of  H  —  xl  =  DAD  —  xl. 

1.  Permute  the  rows  and  columns  of  A  —  iD~^  (which  has  the  same  inertia  as  H  —  xl) 
and  partition  it  as 

Au-xD^^  A, 2 

A21  A22  ~  xD'2 

so  that  if  I  —  xd~^  is  a  diagonal  entry  of  A\\  —  xD^   ,  then  xd~^  >  2n  +  1,  where  n 
is  the  dimension  of  H . 

2.  Compute  X   —  A22  -  xD^^  —  A2\[A\\  —  xD\^)~^ A\2,   using  Cholesky  to  compute 
{Au-xD-^)-'An- 

3.  Compute  inertia(A')  =  [neg,  :ero,pos)  using  a  stable  pivoting  scheme  such  as  in  [4j. 

4.  The  inertia  of  H  —  xl  is  (neg  +  dim(.4ii),  z,p}. 

The  proof  of  correctness  requires  the  following 

Lemma  5.1  Let  H  =  D'A'D'  be  positive  definite,  and  let  Hx  =  b  be  solved  by  Cholesky 
to  get  an  approximate  solution  x.  We  do  not  assume  A'  has  a  unit  diagonal.  Let  £  be  the 
machine  precision,  and  assume  no  overflow  nor  underflow  occurs.   Then  to  first  order  in  e, 

||X  -  ilb   <  0(£M'||2  ■  ||/l'-l||M|Z?'-'||M|6||2) 

Proof.  We  begin  by  defining  some  convenient  notation.  Let  H  be  defined  by  Hi^  = 
(HttHjj)^!'^.  Let  \E\  denote  the  matrix  of  absolute  values  of  entries  of  £^,  and  let  inequalities 
like  X  <  Y  between  matrices  be  interpreted  componentwise.  Then  Lemma  4.14  of  the  last 
section  says  that  if  L  is  the  computed  Cholesky  factor  of  H ,  then  LL^  =  H  +  E  where 
l^"!  <  (n  +  5)en.  *Note  also  by  Cauchy-Schwartz  that  \L\  ■  |I^|  <  H. 

We  begin  by  proving  that  D'(A'+F)D'x  =  b  where  ||F||2  =  0(£-)||A'||2.  In  solving  Ly  =  b 
with  forward  substitution,  we  actually  get  {L  +  6L\  )y  =  b,  where  |^ii,,j|  <  n£\L,j\  [17].  In 
solving  L^x  =  y  we  actually  get  (L  +  fiL2)^x  =  y  where  |<5l2,ij|  <  ne\L^j\.  Altogether 

{H  ^-E^  SLiL'^  +  L6LI  +  8L^6LI)x  =  D'{A'  +  F)D'x  =  b 

where 

IIFII2     =  \\D'-'^{E  +  6LiL'^  ■irLdLl^6Lx6Ll)D'-'^\\2 

<  \\D'-'ED'-'\\2  +  \\D'-'\6Lx\  ■  \L'^\D'-'\\2  +  \\D'-'\L\  ■  |^4|Z)'-'||2  +  \\D'-'\6L,\  ■  \8LI 

<  (n  +  b)e\\D'-^HD'-^\\2  +  n£\\D'-^  H  D'-^\\2  +  n£\\D'-^  H  D''^^  +  n^e'^WD'-^  HD'-^h 

<  {3n'^  +  5n+  n^s)£\\A'\\2 
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Thus 


\x-i\\2     =     \\D'-'A 


'-'  ''-'F(.-l'  +  F)-^D'-^b\\2  <  \\D 


'-^"M|F||2-||.4'-'||^ 


<     0{e)\\D'-%-\\A'\\2-\\A'-'\\l 


to  first  order  in  £. 


I 


Theorem  5.2  Lets  be  the  machine  precision  in  which  Algorithm  5.1  is  carried  out,  where 
we  assume  neither  overflow  nor  underflow  occur.  Then  Algorithm  5.1  computes  the  exact 
inertia  of  D(A  +  8.A)D  -  xl.  where  \\SA\\2  =  0(e).  Thus.  Algorithm  5.1  can  he  used  in 
a  bisection  algorithm  to  find  all  the  eigenvalues  of  H  to  the  accuracy  of  Theorem  2.3  or 
Proposition  2.5. 

Proof.    A'  is  defined  so  tlaat 


/     (.4„-xZ)-2)-i.4i2 


■  .4„  -  xDi'^            Ai2 

A2\                 ^22  —  XD2 

= 

/ 

A2i{Aii-xDi^)-' 

0  ' 
/ 

■  .\ii-xD-^ 
0 

0 
A' 

so  that  the  inertia  of  H  —  il  equals 


inertia(/l  -  xD   ^)  =  inertia(A')  +  inertia(.4n  -  xDj  ^)  =  inertiaiA")  +  (dim(Aii),0,0) 

by  Sylvester's  Theorem  and  the  fact  that  An  -  xDi^  is  negative  definite.  The  algorithm 
in  [4]  will  compute  the  exact  inertia  of  X  +  6X,  where  ||<5A'||2  =  0(£)||A'||2.  Thus  if 
we  show  ||A'||2  =  0(1)  and  that  A'  can  be  computed  with  error  0(e),  we  will  be  done. 
By  construction  the  diagonal  entries  of  .4ii  —  xD^^  are  less  than  or  •  -iial  to  — 2n,  and 
the  ofFdiagonzd  entries  of  all  of  A  -  xD~^  are  bounded  by  1  in  absolute  value.  Write 
-{An  -  xDi'^)  =  Hi  =  DiAiDi  where  Au,  =  1-  Then  ||Ai||2  <  3/2,  ||Af^||2  <  2  and 
||£)^^||2  <  l/(2n).  By  Lemma  5.1  the  error  in  computing  (An  -  xD^^)~^Ai2  is  bounded 
by  0(e).  Also,  by  construction  ||^i2||2  =  ll^2i||2  <  n,  ||(Aii  -  xZ)f^)~^||2  <  1/n  and 
11^22  -  xD2'^\\2  <  3n  +  1,  so  ||A:||2  <  3u  +  1  +  n'^/n  =  4n  +  1  =  0(1)  as  desired.       I 


Algorithm  5.2  Inverse  iteration  for  computing  the  eigenvector  x  of  a  symmetric  positive 
definite  matrix  H  =  DAD  corresponding  to  eigenvalue  z.  tol  is  a  user-specified  stopping 
criterion. 

1.  We  assume  the  eigenvalue  z  has  been  computed  accurately,  for  example  using  Algo- 
rithm 5.1. 

2.  Choose  a  starting  vector  yo;  set  i  =  0. 

3.  Compute  the  LDL^  factorization  of  P{A-zD~^)P^ ,  where  P  is  the  same  permutation 
as  in  Algorithm  5.1,  step  1. 

4-   Repeat 

1  =  1  +  1 
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Solve  (A  -  :D~^)y,  =  (/,_i  for  y,   using  the  LDL^  factorization  of  step  2. 

y,  =  r  ■  y, 
until  (r  <  tol)) 
5.  X  =  D'^y, 

Theorem  5.3  Suppose  Algorithm  5.2  terminates  u-ith  x  as  the  computed  eigenvector  of 
H  =  DAD.  Then  there  is  a  diagonal  matrix  D  with  £>,,  =  1  +0{tol)  and  a  mati'x  dA  with 
|(<5^||2  =  O(tol),  such  that  Dx  is  the  exact  eigenvector  of  D(A  +  SA)D.  Thus,  tue  error  m 
X  IS  bounded  by  Theorem  2.7,  Corollary  2.9  and  Proposition  2.12. 

The  proof  is  identiced  to  the  proof  of  Theorem  11  in  [2]. 
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6      Upper  Bounds  for  max^  k(.4,„)/k(.4o) 

As  stated  in  sections  3  and  4.  our  claims  about  the  accuracy  to  which  Jacobi  can  solve  the 
eigenproblem  depend  on  the  ratio  mdiXm  ^{Am)/^{Ao)  being  modest.  Here  Hq  =  DqAqDo 
is  the  initial  matrix,  and  Hm  =  DmAmDm  is  the  sequence  produced  by  Jacobi  (Hm+\ 
is  obtained  from  Hm  by  applying  a  single  Jacobi  rotation.  Dm  is  diagonal  and  Am  has 
ones  on  the  diagonal).  The  reason  is  that  the  error  bounds  for  Jacobi  are  proportional  to 
maXn,  K(.4m  ).  and  the  error  bounds  of  section  2  are  proportional  to  k(.4o)- 

In  this  section  we  present  several  results  explaining  why  maXm  ^i-^m  )/'<^(-'lo)  should 
not  be  expected  to  grow  very  much.  Recall  that  convergence  of  Hm  to  diagonal  form  is 
equivalent  to  the  convergence  of  Am  to  the  identity  matrix,  or  of  ^c(/lm)  to  1.  Thus  we 
expect  K(Am)  <  k(4o)  eventually.  The  best  situation  would  be  monotonic  convergence, 
but  this  is  unfortunately  not  always  the  case. 

We  have  not  been  able  to  completely  explain  the  extremely  good  numerical  results  of 
section  7,  that  max^  ^(.4^  )/k(.4o)  never  exceeded  1.82.  and  averaged  1.20  in  over  900000 
experiments.  A  complete  theoreticed  explanation  of  this  remains  an  open  question. 

We  will  only  speak  in  terms  of  two-sided  Jacobi  in  this  section.  This  is  no  loss  of 
generality  because  in  exact  arithmetic  one-sided  Jacobi  on  G  is  equivalent  to  two-sided 
Jacobi  on  G^G . 

Our  first  result  will  show  that  K[Am)lfi[Ao)  cannot  be  too  large  if  Am  is  obtained 
from  .4o  by  a  sequence  of  Jacobi  rotations  in  pairwise  disjoint  rows  and  columns.  The 
second  result  give  a  cheaply  computable  guaranteed  upper  bound  on  majCm  «(>im  )/'^(--lo) 
in  terms  of  the  Hadamard  measure  of  .4o.  This  bound  is  generally  quite  pessimistic  unless 
the  dimension  of  .4  is  modest  and  «;(.4o)  is  small,  at  most  a  few  hundred.  The  third  and 
fourth  results  will  be  for  accelerated  one-sided  Jacobi  (Algorithm  4.4).  The  third  result 
shows  that  the  wider  the  range  of  numbers  on  the  diagonal  of  H ,  the  smaller  k[Ai  )  for  that 
algorithm.  This  in  turn  makes  it  converge  faster.  In  other  words,  the  morp  the  guaranteed 
accuracy  of  the  algorithm  exceeds  that  of  QR  (or  any  tridiagonalization  l..i>ed  algorithm), 
the  faster  it  converges.  The  fourth  rather  surprising  result  that  k(A\)  la  bounded  by  a 
constant  depending  only  on  the  dimension  n,  not  on  Ao-  These  last  two  results  lead  us  to 
recommend  accelerated  one-sided  Jacobi  as  the  algorithm  of  choice  (unless  it  is  important 
to  get  small  eigenvector  components  to  high  accuracy;  see  the  discussion  in  subsection  4.3). 

Proposition  6.1  Let  Hq  be  n  by  n.  Let  Hm  be  obtained  from  Hq  by  applying  m  Jacobi 
rotations  in  pairwise  nonoverlapping  rows  and  columns  (this  means  m  <  n/2).  Write 
Hm  =  DmAmDm  OS  before.   Then 

^7^  <min(^(/lo),2n)  (6.2) 

k{Ao) 

Also 

^^^^<min(K(>lo),8)  (6.3) 

Furthermore,  the  spectrum  of  Am  is  independent  of  Dq,  even  though  the  entries  of  Am 
depend  on  Dq.  More  precisely,  the  spectrum  of  Am  coincides  with  the  spectrum  of  the  pencil 
Aq-XA'q,  where  A'q  coincides  with  Ao  on  every  rotated  element  and  is  the  identity  otherwise. 
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Proof.  We  begin  by  deriving  a  matrix  pencil  depending  only  on  .4o  whose  eigenvalues 
are  the  same  as  Am-  This  will  prove  that  the  eigenvalues  of  .4^  depend  only  on  .4o.  We 
assume  without  loss  of  generality  that  the  m  Jacobi  rotations  are  in  rows  and  columns  (1.2). 

(3,4) (2m  -  1.2m).   This  lets  us  write  J^ HqJ  =  H-m  where  J  is  block  diagonal  with 

the  2  by  2  Jacobi  rotations  (and  possibly  ones)  on  its  diagonal.  Rewrite  this  as 

.4„  =  {Dm'J^Do)Ao{DoJDm' )  =  Z^AoZ 

where  Z  has  the  same  block  diagonal  structure  as  J.  Let  ^q  be  a  block  diagonal  matri.x 
with  the  same  block  structure  as  Z  and  J.  where  .4g  is  identical  to  .4o  within  its  2  by  2 
blocks,  and  has  ones  on  its  diagonal  when  J  does.  Since  Hra.\2  =  HmM  —  ■  ■  •  =  0.  also 
Am.\2  =  -4771,23  =  •  •  •  =  0.  Thus  .4^  has  2  by  2  identity  matrices  on  its  diagonal  matching  the 
block  structure  of  Z.  J  and  .4o.  Thus  .4m  =  Z^ AqZ  implies  Z~^ Z~^  =  .4o.  Therefore  the 
eigenvalues  of  .4m  =  Z-^.4o2' are  identical  to  those  of  the  pencil  Aq  -  AZ~'^Z~^  =  ,4o-A^g. 
Now  we  apply  the  minimax  theorem  to  bound  ^miniAm)  below  by 

Amin(-4m)  -  mm    ^  >  =^— —  =  — — 6.4 

li^o  x^  AqI       md^j-^ox^  AqX       1 +  maxi<fc<m  |/lo,2fc-i,2fc| 

We  may  bound  1  +  maxi<fc<m  |-4o.2fc-i,2fc|  from  above  by  both  AmaxC^lo)  and  2.  yielding 

\         C   1      \  -^  Aniin(  Ao)  , 

m\n{2,\m^{AQ)) 

Now  we  bound  Aniax('4m)  from  above.  First  by  the  minimax  theorem  we  may  write 

x'^Aqx        imXi^ox'^Aox        XmnxiAo) 
Aniax(^m)  =  max    J  <  —. 5=——  < — — - 

which  when  combined  with  (6.5)  yields 

proving  half  of  (6.2).    For  the  other  half  note  that  1   <  An,ax(^i)  <  'J  for  all  i,  so  that 
Ainax(^m)/Amax(4o)  <  "■  Now  combine  this  with  (6.5). 

Now  we  show  AmAxC^t+i)  <  4AmAx(^t)i  which  when  combined  with  (6.5)  yields  (6.3).  It 
suffices  to  show  Aniax(^i)  ^  4Aniax(^o)-  Write 


An  = 


An     Ai2 
A21     A22 


where  An  is  2  by  2.  Then  by  the  minimax  theorem  there  exists  a  conformally  partitioned 
unit  vector  x-^  =  [xf  ,12"]  where 


rT"  J..-,.  1  o.,r  j.._.  I  _r 


if  Auxi  4-  2x{  A12X2  +  xj  ^223-2 

'^maxi'^l)  —  T-  1  .       T 

Xf  AiiXi  4-  X^X2 
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Write  xfxi  =  <^  {0  <  (  <  I),  •r2"^2  =  1  -  C-  Jr--lii-ri  =  '^iC  and  x'^ A22^2  =  ■^2(1  -  O-  so 
that 


'^niax(--il)       - 


nC  +  2 jf  .412X2  +  ^2(1  -  C 


The  maximum  of  this  last  expression  over  all  0  <  C  <  1  is 

2  +  2r2  <  2  +  2A^^(.A22)  <  4A^^^(Ao) 

I 

Our  second  bound  is  based  on  the  Hadamard  measure  of  a  symmetric  positive  definite 
matrix  H: 


Proposition  6.6    The  Hadamard  measure  "HiH)  has  the  following  properties: 

1.  H{H)<\  and  n{H)  =  1  if  and  only  if  H  is  diagonal. 

2.  H(H)  =  'H(DHD)  for  any  nonsingular  diagonal  D. 

3.  Let  H  -  DAD  with  D  diagonal  and  A  with  unit  diagonal.    Then 

,      ,  ..^  W^)  _  detU) 
e  e 

where  e  =  exp(l). 

4.  Let  H'  be  obtained  from  H  by  applying  a  Jacobi  rotation  (in  exact  arithmetic)  in  rows 
and  columns  i  and  j.   Then 

n{n')  =  y^>H{n) 

5.  Let  Hq,  . . . ,  Hm^ . . .  be  a  sequence  of  symmetric  positive  definite  matrices  obtained  from 
Jacobi's  method  in  exact  arithmetic.  Let  Hm  =  DmAmDm  with  Dm  diagonal  and  Am 
with  unit  diagonal.   Then 

ne  n-e 

maxK(Am)  < 


det(>lo)       niHo) 

Proof. 

1.  Write  the  Cholesky  decomposition  H  =  LL^ .  Then 

1=1   k=\  t=l 

2.  det(I)^)  factors  out  of  the  numerator  and  denominator  of  HiDHD). 
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3.  From  (2.)  H{H)  =  Hi  A)  =  det(.4).  so  it  suffices  to  show  X.^^JA)  >  deX{A)/e.  Let 
0  <  Ai  <  ■  •  •  <  An  be  the  eigenvalues  of  .-i.  Since  Ai  =  det(.4)/  Y[?=2  •^'-  '■'•'^  ^leed  to  show 
nr=2  ^'  <  e.  Now  ^"=^2-^,  <  tr(.4)  =  n.  Since  ab  >  {a  +  x)(b-  x)  for  all  a  >  6  >  r  >  0,  we 
see  nr=2  ^'  'S  greatest  when  all  A,  =  n/(n-l).  in  which  case  nr=2  ^'  -  (("  -  1 )/")""'  <  e. 

4.  From  Proposition  6.1  we  have 

■H(H')  =  det(.-l.4'-')  =  det(.4)/(l  -  A^)  =  n{H)/{l-  A^) 

where  .4'  =  /  except  for  .4'_,  =  .4^,  -  .4,_,. 

5.  This  is  directly  implied  by  (3.)  and  (4.).      I 

Thus,  Part  5  of  this  proposition  gives  us  a  guaranteed  upper  bound  on  max^  «(.4Tr,)  at 
a  cost  of  about  n'^/6  flops,  compared  to  2n^  flops  per  Jacobi  sweep  {4n^  if  accumulating 
eigenvectors).  If  we  use  the  algorithm  in  subsection  4.3,  where  we  must  do  Cholesky  anyway, 
this  upper  bound  comes  nearly  for  free. 

In  section  7  we  will  present  numerical  results  on  the  utility  of  this  upper  bound.  Basi- 
cally, it  is  only  useful  as  long  as  k(-4o)  is  quite  small  and  Aq  has  low  dimension;  otherwise 
it  is  much  too  large  to  be  useful. 

Our  third  and  fourth  bounds  are  for  accelerated  one-sided  Jacobi  (Algorithm  4.4).  Recall 
that  this  algorithm  begins  by  doing  Cholesky  with  complete  pivoting  on  .^o  to  get  PHqP^  = 
LL^ ,  where  P  is  a  permutation  matrix.  Then  it  does  one-sided  Jacobi  on  L,  which  is 
equivalent  (in  exact  arithmetic)  to  two-sided  Jacobi  ><n.  L^ L.  Therefore,  Algorithm  4.4 
essentially  starts  with  L^ L  =  Hi  =  DiAiDi.  The  trauaition  from  Hq  to  Hi  is,  in  fact,  one 
step  of  the  symmetric  LR  algorithm  which  usually  has  some  non-trivial  diagonalizing  effect 
(the  pivoting  cares  for  the  proper  ordering). 

Our  third  result,  which  we  state  rather  informally,  is  that  the  larger  the  range  of  numbers 
on  the  diagonad  D^  of  H,  the  smaller  is  k(.4i).  This  effect  was  also  observed  in  [16].  We 
argue  as  follows.  Let  L  =  DLa  be  the  factor  obtained  from  complete  pivoting.  Here,  La 
has  rows  of  unit  norm.  Since  Algorithm  4.4  does  one-sided  Jacobi  on  Z,  its  performance 
depends  on  the  condition  number  of  DLaD',  where  D'  is  chosen  diagonal  to  make  the 
columns  of  DLaD'  unit  vectors.  From  van  der  Sluis's  theorem  [13]  we  know  the  condition 
number  of  DLaD'  can  be  at  most  n  times  DLaD~^ ,  so  it  suffices  to  examine  k{DLaD~^). 
The  effect  of  complete  pivoting  is  essentially  to  reorder  D  so  that  D„  >  Z?,+i,,+i,  and  to 
keep  LA,ii  as  large  as  possible.  Now  (DLaD'^)^,  =  La.h  is  unchanged,  and  the  subdiagonal 
entry  {DLaD~^),j  =  L^,,jZ?„Z)j~'  is  multiplied  by  the  factor  Dt,D~^  which  is  between  0 
and  1.  The  more  Djj  exceeds  Z)„,  the  smsdler  this  factor,  and  the  more  nearly  diagonal 
DLaD~^  becomes.  Since  complete  pivoting  tries  to  keep  the  diagonal  of  La  large,  this 
improves  the  condition  number.  The  acceleration  from  this  effect  will  be  quite  pronounced 
in  the  numerical  examples  of  section  7. 

Our  fourth  result  shows  that  surprisingly  maXm>i  ii{Am)  is  bounded  independent  of  Hq: 

Proposition  6.7  Let  PHqP^  =  LL^  be  the  Cholesky  decomposition  of  the  n  by  n  matrix 
Hq  obtained  with  complete  pivoting.  Let  Hx  =  L^ L  =  DiAiDi.  Let  Hm  =  DmAmDm, 
m  >  1,  6e  obtained  from  two-sided  Jacobi  applied  to  Hi.   Then 

1.  n(Hi)>n(Ho). 

2.  'H(Hi)  >  1/n!.   This  bound  is  attainable. 
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3.   m'ax„>iK(.4^)  <  e-n/HiHi)  <  e  ■  n  ■  n\ 

Proof. 

1.  Since  det{Hi)  -  det{Ho),  it  suffices  to  show  H.  ^o.n  >  Hi  ^i.a-  Assume  without 
loss  of  generality  that  P  =  I.  Then  Ho,n  =  T,'k=i  ^h  ^^'^  ^i."  =  ELi-^L-  Complete 
pivoting  is  equivalent  to  the  fact  that  Lf,  >  Yli=,  ^%  ^°^  ^^^  J  -  '■  ^^"^  ^i^''  ^'^  P^ove 
nr=i  Yi\=i  L^k  ^  nr=i  Z?=:  -^i-  ^^'^  systematically  use  the  fact  that  ab  >  [a  +  xj[^b  -  x) 
for  a  >  6  >  J  >  0.  We  illustrate  the  general  procedure  in  the  case  of  n  =  3: 

—       (-^11  +  -^21  "t"  -^31 )( ■'^22)' -^32  +  ■''33) 
>       (Z.)!  +  L21  +  L3j)(Z/22  +  •^32)(-^33) 

2.  We  have 

det(i)2         nr=i  ^?, 


n{H,)  = 


"  £2  "11 


To  see  that  this  bound  is  attainable,  let  H  =  LL^  where  Z„  =  ^''~^'/^  and  L,-,  =  (1  - 

|iji/2^(>-i)/2    >y;Q^  let  /i  >  0  become  small. 

3.  The  result  follows  from  part  2  and  Proposition  6.6,  part  5.       I 

The  example  in  part  2  of  the  Proposition  for  which  the  Hadamard  bound  is  attainable 
unfortunately  has  the  property  that  the  resulting  upper  bound  in  part  3  is  a  gross  overesti- 
mate. While  the  upper  bound  grows  as  e-n-n!,  k{A\)  only  grows  like  r?-''-'  However,  «(Ao) 
grows  like  ^""Z^,  which  can  be  arbitrarily  larger  than  the  bound  in  paii  3.  The  choice 
/x  =  .5  provides  an  example  where  the  upper  bound  in  part  3  can  arbitrarily  exceed  both 
k(^o)  and  maxm>i  K{Am)  for  large  n. 

Nonetheless,  in  numerical  experiments  the  upper  bound  e  ■n/'H(Hi)  on  maxm>i  «(/!„) 
never  exceeded  40.  We  also  always  observed  that  k(Ai)  <  k{Ao)  in  all  rases,  although  we 
have  not  been  able  to  prove  it  in  general. 
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7      Numerical  Experiments 

In  this  section  we  present  the  results  of  numerical  experiments.  Briefly,  we  tested  every 
error  bound  of  every  algorithm  presented  in  this  paper,  and  verified  that  they  held  in  all 
examples.  In  fact,  the  performance  is  better  than  we  were  able  to  explain  theoretically, 
both  because  we  could  observe  little  or  no  growth  in  actual  errors  for  increasing  dimension, 
and  because  of  the  surprisingly  small  values  attained  by  max,„  k(  Am)/K(/lo)  (see  section 
6). 

These  tests  were  performed  using  FORTRAN  on  a  SUN  4/260.  The  arithmetic  was 
IEEE  standard  double  precision  (1],  with  a  machine  precision  of  £■  =  2~^^  ~  10~'^  and 
over/underflow  threshold  10*''°*. 

There  were  essentially  four  algorithms  tested:  two-sided  Jacobi  (Algorithm  3.1),  one- 
sided Jacobi  (Algorithms  4.2  and  4.1).  accelerated  one-sided  Jacobi  (Algorithms  4.4  and 
4.1).  and  bisection/inverse  iteration  (Algorithms  5.1  and  5.2).  All  were  used  with  the 
stopping  criterion  tol  =  lO"'**. 

Since  we  claim  these  algorithms  are  more  accurate  than  any  other,  we  tested  their 
accuracy  as  follows.  We  considered  only  symmetric  positive  definite  eigenproblems,  and 
solved  every  one  using  every  algorithm.  The  different  answers  were  compared  to  see  if 
they  agreed  to  the  predicted  accuracy  (which  they  did).  They  were  also  compared  to 
the  EISPACK  routines  tred2/tqr2  [14],  which  implement  tridiagonalization  followed  by  QR 
iteration.  Small  eigenvalues  computed  by  EISPACK  were  often  negative,  indicating  total 
loss  of  relative  accuracy. 

Then  rest  of  this  section  is  organized  as  follows:  Subsection  7.1  discusses  test  matrix 
generation.  Subsection  7.2  discusses  the  accuracy  of  the  computed  eigenvalues.  Subsection 
7.3  discusses  the  accuracy  of  the  computed  eigenvectors.  Subsection  7.4  discusses  the  the 
growth  of  maxm  k(/1^)/k(  Aq).  Subsection  7.5  discusses  convergence  rates;  here  the  speed 
advantage  of  accelerated  one-sided  Jacobi  will  be  apparent. 

7.1      Test  Matrix  Generation 

We  generated  several  categories  of  random  test  matrices  according  to  three  parameters: 
the  dimension  n,  ka,  and  kq.  First  we  describe  the  algorithm  used  to  generate  a  random 
matrix  from  these  parameters,  and  then  the  sets  of  parameters  used. 

We  tested  matrices  of  dimension  n  =  4,  8,  16  and  50.  Since  testing  involved  solving  an 
n  by  n  eigenproblem  after  each  Jacobi  rotation  (to  evaluate  K(Am))  and  there  are  0{n^) 
Jacobi  rotations  required  for  convergence,  testing  costs  0{n^)  operations  per  matrix. 

Given  k^i,  we  generated  a  random  symmetric  positive  definite  matrix  with  unit  diagonal 
and  approximate  condition  number  ka  as  follows.  We  began  by  generating  a  diagonal  matrix 
T  with  diagonal  entries  in  a  geometric  series  from  1  down  to  1/ka-  Then  we  generated  an 
orthogonal  matrix  U  uniformly  distributed  with  respect  to  Haar  measure  [15],  and  formed 
UTU^ .  Finally,  we  computed  another  diagonal  matrix  K  so  that  Aq  -  KUTU^K  had 
unit  diagonal.  This  last  transformation  can  decrease  the  condition  number  of  UTU'^ ,  but 
usually  not  by  much.  For  4  by  4  matrices,  it  decreased  it  by  as  much  as  a  factor  of  500, 
for  8  by  8  matrices  by  a  factor  of  20,  for  16  by  16  matrices  by  a  factor  of  5  and  for  50  by 
50  matrices  by  a  factor  of  1.5.  (This  decreasing  variability  is  at  least  partly  due  to  the  fact 
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that  we  ran  fewer  tests  on  the  larger  matrices.)  For  a  more  complete  discussion  of  the  test 
matrix  generation  software,  see  [8]. 

Given  kd«  we  generated  a  random  diagonal  matrix  Do  with  diagonal  entries  whose 
logarithms  were  uniformly  distributed  between  0  and  logK^j.  This  means  the  diagonal 
entries  themselves  were  distributed  from  1  to  Kp-  The  uniform  distribution  of  the  logarithm 
essentially  means  every  decade  i-  •  qually  likely,  and  so  generates  matrices  Dq  w-ith  entries 
of  widely  varying  magnitudes. 

The  resulting  random  matrix  was  then  Hq  =  DqAoDq. 

We  generated  random  matrices  with  5  possible  different  vcilues  of  K.4:  10,  lO'^.  IC*,  10® 
and  10^2  6  possible  different  values  of  Kp:  10^  10^°.  lO^o,  lO^o,  iQ^o  and  10i°°.  and  4 
different  dimensions  n  =  4.  8.  16  and  50.  This  makes  a  total  of  5  x  6  x  4  =  120  different 
classes  of  matrices.  In  each  class  of  dimension  n  =  4  matrices,  we  generated  100  random 
matrices,  in  each  class  of  n  =  8.  we  generated  50  random  matrices,  in  each  class  of  n  =  16, 
we  generated  10  random  matrices,  and  in  each  class  of  n  =  50,  we  generated  one  random 
matrix.  This  makes  a  total  of  4830  different  test  matrices. 

The  matrices  had  in  some  cases  eigenvcdues  ranging  over  200  orders  of  magnitude  (when 
KD  =  10^°°).  The  relative  gaps  relgapx  ranged  from  .028  to  2  •  lO''^ 

7.2      Accuracy  of  the  Computed  Eigenvalues 

There  are  two  accuracy  bounds  for  eigenvalues  from  section  2  which  we  tested.  The  first 
one  is  based  on  Theorem  2.3  (or  Theorem  2.17  together  with  Theorem  4.15),  which  says 
that  if  A(  and  A"  are  approximations  of  A,  computed  by  two  of  our  algorithms,  then 

|A'  -  A"| 


k(.4o)A; 


should  be  0{tol),  where  tol  =  lO"^**  is  our  stopping  criterion.  For  two-sided  Jacobi  and 
one-sided  Jacobi,  Qi  never  exceeded  2  •  10"^^.  For  two-sided  Jacobi  and  accelerated  one- 
sided Jacobi,  Q\  also  never  exceeded  2  •  10"^^.  Every  matrix  had  an  eigenvalue  for  which 
Qi  exceeded  4  •  10"^*,  showing  that  the  bound  of  Theorem  2.3  is  attainable,  as  predicted 
by  Proposition  2.13. 

In  the  case  of  bisection,  we  did  not  run  a  bisection  algorithm  tn  convergence  for  each 
eigenvalue,  but  rather  took  the  eigenvalues  X[  computed  by  two-sided  .Jacobi,  made  intervals 
[(1  -  tol  ■  k(Ao))A(,  (1  +  tol  ■  k(Ao))X[]  from  each  one,  and  used  bisection  to  verify  that 
each  interval  contained  one  eigenvalue  (overlapping  intervals  were  merged  and  the  counting 
modified  in  the  obvious  way).  All  intervals  successfully  passed  this  test. 

The  second  accuracy  bound  is  from  Proposition  2.5  (or  Proposition  2.19  together  with 
Theorem  4.15)  which  predicts  that 

I  A'  -  Al'l 


Q2=    nr.    ..,12 


li^oi'.l 


2 


should  be  0{tol).  Here  v,  is  the  unit  eigenvector  computed  by  two-sided  Jacobi.  For  two- 
sided  Jacobi  and  one-sided  Jacobi,  Q2  never  exceeded  2  •  lO"^"*.  For  two-sided  Jacobi  and 
accelerated  one-sided  Jacobi,  Q2  never  exceeded  9  •  10"^^.  Every  matrix  had  an  eigenvalue 
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for  which  Q2  exceeded  5  •  10  ^^ .  showing  that  the  bound  of  Proposition  2.5  is  attainable, 
as  it  predicts. 

In  the  case  of  bisection,  we  again  made  intervals  [A;  -  tol  ■  \\Dov^\\l,\[  +  tol  ■  \\Dqi\\\\] 
from  each  eigenvalue  A'  and  verified  that  each  interval  contained  the  proper  number  of 
eigenvalues. 

Finally,  we  verified  a  slightly  weakened  version  of  Proposition  2.10.  that 

'^min(-lo)  -tol<-^<  X^^^iAo)  +  tol 

n, 

for  the  eigenvalues  A^  computed  by  two-sided  Jacobi.  Here  h,  is  the  i'-th  imallest  diagonal 
entry  of  Hq-  .Adding  and  subtracting  tol  to  the  upper  and  lower  bounds  takes  into  account 
the  errors  in  computing  X[. 

7.3     Accuracy  of  the  Computed  Eigenvectors 

There  is  one  bound  on  the  magnitude  of  the  components  of  the  eigenvectors,  and  two 
accuracy  bounds,  one  for  the  norm  error  and  one  for  the  componentwise  error. 

We  begin  with  a  few  details  about  our  implementation  of  inverse  iteration.  We  used  the 
eigenvalues  computed  by  two-sided  Jacobi,  and  the  vector  of  all  ones  as  a  starting  vector. 
Convergence  always  occurred  after  just  one  iteration. 

The  componentwise  bound  on  the  magnitude  of  the  eigenvectors  is  based  on  Proposition 
2.11,  which  says  that  the  components  of  the  normalized  eigenvector  i;,  should  be  bounded 

by 

/  A   \^^^       /  \    \l/2 
■      \vdJ)\<vdj)  =  iK{Ao)f^'-nnn{i^^]        ,  ^^j      ) 

This  was  verified  for  the  eigenvectors  computed  by  all  four  algorithms.  We  note  that  since 
this  bound  is  proportional  to  k(Ao)'^^'^,  it  becomes  weaker  as  k{Ao)  becomes  larger,  and 
indeed  becomes  vacuous  for  matrices  with  k{Ao)  large  and  eigenvalues  in  a  narrow  range. 
The  norm  error  bounds  are  based  on  Theorem  2.7  (or  Theorem  2.21  together  with 
Theorem  4.15),  which  predicts  that  if  v'^  and  v"  are  approximations  of  the  unit  eigenvector 
Vi  computed  by  two  of  our  algorithms,  then 

„    _  llv.'-r."||, 

Q3  = 


{K{Ao)/relgapx,)  +  1 

should  be  O(tol).  (We  add  the  1  in  the  denominator  because  a  single  roundoff  error  in 
the  largest  entry  can  cause  a  norm  error  of  e;  see  Theorem  3.11  or  Theorem  4.9.)  For 
two-sided  Jacobi  and  one-sided  Jacobi,  Q3  never  exceeded  3  •  10"^^.  For  two-sided  Jacobi 
and  accelerated  one-sided  Jacobi,  Q3  also  never  exceeded  2  •  lO"^**.  For  two-sided  Jacobi 
and  inverse  iteration,  Q3  never  exceeded  8  •  IG"^"*.  Every  matrix  had  an  eigenvector  for 
which  Q3  exceeded  10"^*  for  every  pair  of  algorithms  compared,  showing  that  the  bound 
of  Theorem  2.7  is  nearly  attainable,  as  predicted  by  Proposition  2.14. 

The  second  accuracy  bound  is  based  on  Proposition  2.12  (or  Proposition  2.26  and  The- 
orem 4.15),  which  predicts 

g^  _  KU)  -  <(j)tmin(re/gapA.,2-i/^) 
k{Ao)  ■  u,(j) 
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Table  1: 

Hadamard  upper 

bound  Qe 

on  maxm  «(^m)/«(--lo) 

n 

^A 

10 

10^ 

lO'* 

10»                   10^2 

4 

5.8 

13 

590 

6.3  ■lO*'          6.1 -lO^*^ 

8 

21 

410 

1.1-10" 

9.1-10''                      oo 

16 

200 

2.7  •  10^ 

1.8-  10^^ 

oc                      cc 

50 

6.4 

•  10^ 

8.0  •  lO^*^ 

cc 

oc                       oc 

should  be  O(tol).  For  two-sided  Jacobi  and  one-sided  Jacobi.  Q4  never  exceeded  3  •  10"'". 
For  two-sided  Jacobi  and  inverse  iteration,  Q4  never  exceeded  3-10"'^.  For  two-sided  Jacobi 
and  accelerated  one-sided  Jacobi.  Q4  was  as  large  as  .02,  which  is  consistent  with  the  fact 
that  accelerated  one-sided  Jacobi  computes  the  eigenvectors  as  left  singular  vectors  of  L, 
for  which  we  only  have  a  normwise  error  bound  (Theorem  4.9).  For  the  other  algorithm 
Q4  was  only  10"^°  for  matrices  with  k{Ao)  =  10'^;  this  reflects  the  factor  k(/1o)^''^  in  the 
denominator  of  Q^,  a  weakness  of  Proposition  2.11.  In  other  words,  the  componentwise 
error  bounds  are  generally  only  interesting  for  small  to  medium  k,(Ao). 

7.4      Growth  of  maXr7i-«(-4m)/«(-4o) 


In  computing 


Q5  =  rmxK(Am)/K(Ao) 


we  note  that  a  single  computation  requiring  M  Jacobi  rotations  supplied  us  not  just  with 
one  value  of  Q5  but  rather  A/  -  1:  Since  every  ^4,  can  be  thought  of  as  starting  a  new 
eigenvalue  computation,  we  may  also  measure  maxm>,  K{Am)/'<-(^i)  for  all  i  <  M.  Thus, 
all  told,  our  4830  different  matrices  represent  over  900000  data  points  of  Q5. 

The  largest  value  of  Q^  encountered  was  1.82.  This  was  for  an  8  by  8  matrix  with 
k{Aq)  =  1.4  •  10^2,  and  eigenvalues  ranging  over  133  orders  of  magnitude.  141  Jacobi 
rotations  (a  little  over  5  sweeps)  were  required  for  convergence,  plus  28  more  steps  (one 
more  sweep)  where  no  work  is  done  to  recognize  convergence.  In  Figure  1,  a  plot  is  shown 
of  K(Ai)  —  1  versus  i.  We  plot  k{A,)  -  1  instead  of  k(A,)  in  order  to  see  the  quadratic 
convergence  of  k{A,)  to  1.  The  graph  appears  nearly  monotonic,  except  for  a  slight  rise 
near  i  =  20.  This  is  seen  more  clearly  in  Figure  2,  which  plots  maXm>t  «(^tn)/'^('4,)  versus 
i.  Here  the  maximal  nonmonotonicity  of  the  curve  near  i  =  20  is  apparent. 

Now  we  consider  the  Hadamard  based  upper  bound  on  Qs  from  Proposition  6.6: 


Qs  <  C?6  = 


e  ■  n 


niHo)-K(Ao) 


Table  1  gives  the  maximum  values  of  this  upper  bound  for  different  values  of  dimension 
n  and  ka  *  k(Ao).  Recall  that  the  true  value  of  Q5  never  exceeds  1.82.  As  Proposition 
6.6  suggests,  this  upper  bound  should  not  depend  on  Dq  and  indeed  the  values  observed 
depended  very  little  on  Dq. 

As  can  be  seen,  the  Hadamard  based  bound  is  of  little  use  except  for  very  small  matrices 
of  modest  k{Ao).  00  means  the  value  overflowed. 
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Figure  1:  k(.4,)  -  1  versus  i 
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Figure  2:  ma.Xm>,  h{A„.)/k{A,)  versus  i 
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Now  we  consider  accelerated  one-sided  Jacobi.  Let  us  rei .  ne  notation  of  section  6: 
Let  PHqP^  =  LL'^  be  Cholesky  with  complete  pivoting,  and  let  L'^ L  =  Hi  =  DiAiDi. 
As  suggested  in  that  section,  we  expect  both  k(.4i)  to  be  smallf^r  than  k(.4o).  and  the 
Hadamard  based  upper  bound 

^                              s  ■  n 
Qs  <  Qt  =  max  1,- — ) 

on  Qi  to  be  much  smaller  than  the  one  for  two-sided  Jacobi. 

First  of  all  k{Ai)/k{Ao)  never  e.xceeded  .6.  In  fact.  k(Ai)  never  exceeded  4O  for  any 
matrix.  This  is  quite  remarkable.  This  means  that  all  essential  rounding  errors  occurred 
during  the  initial  Cholesky  decomposition.  Finally,  the  Hadamard  upper  bound  Q-  on  Q.5 
never  e.xceeded  29.  This  leads  us  to  strongly  recommend  using  the  Hadamard  measure  in 
conjunction  with  accelerated  one-sided  Jacobi  to  get  error  bounds. 

7.5      Convergence  Rates 

We  begin  with  a  few  details  of  how  we  counted  the  number  of  Jacobi  rotations  required  for 
convergence.  In  all  three  algorithms  (two-sided  Jacobi,  one-sided  Jacobi  and  accelerated 
one-sided  Jacobi),  we  stopped  when  the  last  n{n-l)/2  stopping  tests  \H,j\-{H„Hjj  )^^^  <  tol 
succeeded;  this  means  every  offdiagonal  entry  of  H  satisfies  the  stopping  criterion.  In  the 
case  of  two-sided  Jacobi,  this  means  the  last  71(71  —  l)/2  Jacobi  rotations  involved  almost 
no  work.  For  the  two  one-sided  Jacobis,  however,  evaluating  the  stopping  criterion  costs  3 
inner  products,  so  the  last  n[n  -  \.)l'2  rotation  involve  a  significant  amount  of  work,  even 
if  no  rotations  are  performed.  This  must  be  kept  in  mind  when  comparing  the  number  of 
rotations  for  two-sided  and  one-sided  Jacobi. 

We  used  the  same  standard  cyclic  pivot  sequence  for  all  the  algorithms:  (1,2),  (1.3),  .... 
(l,n),(2,3) (2,7i),(3,4),  ...,  (7i-l,n). 

We  begin  by  comparing  two-sided  Jacobi  and  one-sided  Jacobi.  li  v;act  arifhmetic, 
these  two  algorithms  are  identical.  In  practice,  they  usually  took  the  saiii'  numbei  if  steps, 
although  one-sided  Jacobi  did  vary  from  20%  faster  to  50%  slower  than  two-sided  Jacobi 
on  some  examples.  From  now  on  we  will  only  compare  two-sided  Jacobi  to  accelerated 
one-sided  Jacobi. 

The  most  interesting  phenomenon  was  the  speed  up  experienced  by  accelerated  one- 
sided Jacobi  with* respect  to  two-sided  Jacobi.  In  Table  2  we  present  the  raw  data  on  the 
number  of  sweeps  required  for  convergence. 

There  are  a  number  of  interesting  trends  exhibited  in  this  table.  First,  AOsJ  (accelerated 
one-sided  Jacobi)  never  takes  more  than  6  sweeps  to  converge  for  any  matrix,  wher  -  TsJ 
(two-sided  Jacobi)  takes  up  to  16.5.  In  fact  AOsJ  is  almost  ah\:i\-  faster  than  1-1  (in 
one  example  it  took  5%  longer),  and  can  be  up  to  5  times  faster  (3.0  - v'^eps  vs.  15.2 
sweeps  for  k^  =  10^^,  kd  =  10^*^  and  n  =  50).  Second.  1  he  number  of  sweeps  !'  ftpases 
with  increasing  ka  for  TsJ,  but  not  for  AOsJ.  Third,  the  number  of  sweeps  incri.  -  with 
increasing  dimension  for  both  TsJ  and  AOsJ,  but  much  more  m^Jestl}  for  AOsJ  (from 
2-3  up  to  6)  than  for  TsJ  (from  3-4  up  to  15).  Thus,  the  running  time  for  AOsJ  is  much 
less  dependent  on  the  problem  size  or  sensitivity  (as  measured  by  k^)  than  TsJ.  Fourth, 
the  number  of  sweeps  descreases  as  kq  increases,  both  for  TsJ  and  AOsJ,  but  much  more 
markedly  for  AOsJ  (up  to  a  factor  of  2)  than  for  TsJ  (usually  just  1  sweep). 
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Table  2:  Average 

Number  of  S 

^veeps  for  Two- 

sided  Jacobi  (TsJ) 

and  A 

ccelerat 

ed  One-sided  Jacobi 

(AOsJ) 

ka 

KD 

Dimension  n 

4 

8 

16 

dO 

TsJ 

AOsJ 

TsJ 

AOsJ 

TsJ 

AOsJ 

TsJ 

AOsJ 

10 

10^ 

3.7 

3.0 

4.9 

3.7 

5.7 

4.4 

6.4 

5.0 

lO^o 

3.5 

2.5 

4.6 

3.3 

5.6 

4.1 

6.4 

5.0 

lO^o 

3.1 

2.2 

4.5 

2.8 

5.5 

3.6 

6.0 

4.0 

103O 

3.0 

2.1 

4.6 

2.5 

5.5 

3.4 

6.3 

4.0 

IQSO 

2.8 

1.9 

4.4 

2.3 

5.5 

3.1 

5.8 

4.0 

IQlOO 

2.7 

1.7 

4.5 

2.0 

5.6 

2.6 

5.8 

3.0 

10'^ 

10^ 

3.8 

3.0 

5.2 

3.8 

6.4 

4.5 

7.5 

6.0 

10^° 

3.5 

2.5 

5.1 

3.3 

6.2 

4.1 

7.4 

5.0 

1020 

3.2 

2.2 

4.9 

2.9 

6.2 

3.9 

7.1 

4.0 

103O 

3.0 

2.0 

4.8 

2.6 

5.8 

3.3 

6.8 

4.1 

Iqso 

2.9 

1.9 

4.8 

2.2 

6.1 

3.0 

6.5 

4.0 

IQlOO 

2.8 

1.6 

4.7 

2.0 

6.0 

2.7 

6.8 

3.4 

lO'' 

10^ 

4.0 

2.9 

5.8 

3.6 

7.5 

4.5 

9.2 

6.0 

10^° 

3.7 

2.5 

5.6 

3.3 

7.2 

4.1 

9.3 

5.0 

lO^o 

3.2 

2.2 

5.3 

2.9 

7.2 

3.7 

8.5 

4.9 

103O 

3.1 

2.1 

5.2 

2.6 

6.8 

3.1 

8.2 

4.0 

1050 

2.9 

1.9 

5.2 

2.4 

6.6 

3.0 

8.5 

4.6 

IQlOO 

2.7 

1.7 

4.9 

2.2 

6.9 

2.4 

8.0 

3.9 

10** 

10^ 

3.9 

2.7 

6.4 

3.5 

9.7 

4.1 

13.5 

6.0 

IQio 

3.6 

2.3 

6.3 

3.2 

9.4 

3.8 

12.4 

5.0 

102° 

3.3 

2.1 

5.7 

2.8 

8.9 

3.5 

11.7 

4.7 

1030 

3.1 

2.1 

5.5 

2.6 

8.6 

3.4 

12.0 

4.0 

•10^° 

2.9 

1.9 

5.3 

2.3 

8.5 

3.1 

11.6 

4.0 

IQlOO 

2.9 

1.7 

5.1 

2.0 

8.7 

2.6 

11.6 

4.0 

IQi'^ 

10^ 

3.8 

2.5 

6.8 

3.1 

10.6 

4.0 

16.5 

6.0 

101° 

3.6 

2.2 

6.4 

3.0 

10.3 

3.9 

15.6 

5.0 

1020 

3.4 

2.1 

6.0 

2.7 

9.8 

3.5 

15.3 

5.0 

103O 

3.1 

2.0 

5.8 

2.5 

10.2 

3.3 

15.2 

4.0 

IQSO 

2.9 

1.9 

5.6 

2.3 

9.3 

3.2 

13.7 

3.9 

IQlOO 

2.8 

1.6 

5.2 

2.0 

8.7 

2.7 

15.2 

3.0 
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8      Conclusions 

In  this  paper  we  have  developed  new  perturbation  theory  for  the  eigenvalues  and  eigenvec- 
tors of  symmetric  positive  definite  matrices,  as  well  as  for  eigenvalues  of  symmetric  positive 
definite  pencils.  This  theory  assumes  the  perturbations  are  scaled  analogous  to  the  way 
the  matri.x  is  scaled,  letting  us  derive  much  tighter  bounds  than  in  the  classical  theory.  In 
particular,  we  get  relative  error  bounds  for  th'-  eigenvalues  and  individual  components  of 
the  eigenvectors.  whi(  li  .ire  (nearly)  attainable.  The  bound  for  symmetric  positi%'e  definite 
pencils  may  be  applied  to  matrices  arising  in  finite  element  modeliiiu. 

Second,  we  have  shown  both  through  formal  error  analysis  and  numerical  experiment 
that  Jacobi"s  method  (with  a  modified  stopping  criterion)  computes  the  eigenvalues  and 
eigenvectors  with  these  error  bounds.  We  also  show  that  bisection  and  inverse  iteration 
(applied  to  the  original  matrix)  attain  these  bounds.  In  contrast,  methods  based  on  tridi- 
agonalization  (such  as  QR,  divide  and  conquer,  traditional  bisection,  etc.)  fail  to  attain 
these  bounds. 

We  have  similar  perturbation  theorems  for  the  singular  valuf^  doromposiiion  of  a  general 
matrix  and  the  generalized  singular  values  of  a  pair  of  matrices,  and  similar  error  analyses 
and  numerical  experiments  for  one-sided  Jacobi  applied  to  tlil^  problem.  We  may  also  use 
one-sided  Jacobi  to  solve  the  symmetric  positive  definite  eigenproblem. 

We  have  discussed  an  accelerated  version  of  Jacobi  for  the  «;vmmetric  positive  defi- 
nite eigenproblem,  which  has  the  property  that  the  more  its  accuiacy  exceeds  that  of  QR 
(or  other  conventioncd  algorithms),  the  faster  it  converges.  However,  it  cannot  '-impute 
tiny  components  of  eigenvectors  as  accurately  as  the  other  versions  of  Jacobi,  although  it 
computes  the  eigenvectors  with  the  same  norm  error  bounds.  Unless  getting  the  tiny  eigen- 
vector components  is  important,  we  recommend  this  accelerated  version  of  Jacobi  for  the 
symmetric  positive  definite  eigenproblem. 

The  quantity  max^  K(.4m)/«(,4o)  was  seen  to  be  central  in  the  analysis  of  Jacobi's  accu- 
racy. Numerical  experiments  show  it  to  be  much  smaller  in  practice  than  we  can  explain.  For 
the  accelerated  version  of  Jacobi  we  provide  an  inexpensive  estimator  of  max^  k(  Am)/«(^o) 
which  works  very  well  in  practice.  Explaining  the  excellent  behavior  of  maXm  K{Am  )/k{Ao) 
is  an  important  open  problem. 

The  error  analyses  of  Jacobi  dealt  only  with  the  simplest  implementations.  It  would  be 
worthwhile  to  extfnd  these  analyses  to  cover  various  enhjincements  introduced  by  Rutishauser, 
Veselic,  Hari  and  others.  These  include  delayed  updates  of  the  diagonal  entries  and  an  alter- 
nate formula  for  updating  the  offdiagonal  entries  [12,16],  as  well  as  block  Jacobi  methods. 

In  future  work  we  plan  to  extend  these  results  to  the  symmetric  positive  definite  gener- 
alized eigenproblem,  as  weU  as  indefinite  matrices.  Any  extension  requires  an  appropriate 
perturbation  theory;  therefore  we  do  not  expect  to  be  able  to  extend  the  result  to  jdl  indefi- 
nite matrices,  since  there  is  no  guaranteed  way  to  compute  the  zero  eigenvalues  of  a  singular 
matrix  to  "high  relative  accuracy"  without  computing  them  exactly,  a  feat  requiring  high 
precision  arithmetic.  A  class  of  indefinite  matrices  for  which  a  suitable  perturbation  theory 
exists  are  the  scaled  diagonally  dominant  matrices  [2].  The  perturbation  theory  also  already 
exists  (at  least  for  eigenvalues)  for  the  symmetric  positive  definite  generalized  eigenproblem 
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